EEG Analyses

One of the goals of the RDoC project is to relate measures of working memory across modalities. In the fMRI data, we found differences in load effects based on working memory capacity - typically an inverted U shape relationship between capacity and load effects. Here, we look to investigate whether there are any siimilar relationships in the EEG data, in addition to attempting to replicate some findings in the literature.

ERP data:

During both the cue and probe, we see significant load effects in the N170 components in the average of the O1 and O2 electrodes and in the average of PO7 and PO8 (which is slightly stronger). For this effect, there’s a numerical inverted U shape relationship with capacity. We also see a load effect in the P3 in the parietal electrode average; for this effect, we see a linear increase with capacity. Both of these show a larger amplitude in the low load than high load.

When we look at differences across WM capacity groups, we see a linear increase in activity in P3 in the high load condition.

Spectral data:

Mid occipital cluster shows a positive correlation with span during the encoding and delay period in the beta and low gamma bands, and a negative correlation in the theta through alpha bands during the entire task, but most prominently right at the beginning of the probe period.

Oz shows a strong positive correlation with alpha power during encoding and what looks to be low gamma across the entire task. There also seems to be a negative correlation with beta across the entire task, with particularly strong negative correlations right at the end of the delay period/beginning of probe period.

F avg shows correlations during encoding with the theta band and low gamma/high beta that drift into the beginning of delay, and then a strong negative correlation right at the end of the delay period/beginning of the probe period.

There is also a strong correlation with the fMRI load effect ini the alpha/beta range across all electrodes, especially in encoding.

When we look at all subjects, we see load effects during the cue period in: - alpha, beta: Oz, mid occip, F avg (low > high) - theta: Oz electrode (high > low)

delay period: - alpha, theta: Oz electrode (high > low) - beta: Oz and mid occipital (high > low) - low gamma: mid occipital, Oz (low > high)

probe period; - alpha, beta, theta, low gamma: Oz, mid occip, F avg (low > high)

When we split spectral data into WMC groups, we see:

cue period: - theta: Oz electrode show larger electrodes in high WMC subjects - low gamma: mid occipital (high load and LE) and Oz show med WMC > high WMC (high load and LE) and mid occipital also shows low > high in raw power; trend towards inverted U shaped relationship for LEs in Oz electrode

delay period: - alpha: (non-significant) inverted U shape relationship in Oz - theta: mid occipital cluster shows negative relationship between capacity and LE and high load; more negative power in raw high load data in high WMC > low WMC - low gamma: mid occipital, Oz show low WMC > high WMC and Oz shows med WMC > high WMC in raw power during high load trials mid occipital, Oz electrode show low WMC > high WMC load effects

probe period: - alpha: mid occipital cluster high load and load effect - high WMC more neg raw power than low WMC and more neg LE than med WMC

For the high load, we see a significant relationship between WMC group and theta-beta coupling in the O2 electrode, with a positive linear trend. We also see trends in the Fz electrode for the theta-beta (inverted U shape relationship) and theta-hgamma PAC (negative linear trend). For the load effects, the only significant PAC measure across WMC groups is the theta - beta coupling in the Fz electrode, with an inverted U shaped relationship.

CDA only showed a relationship with the LCD K max measure in the 1 dot condition, but no correlation or obvious inverted U with the DFR load effects or our omnibus span measure. Numerically, it looks like there is an inverted U shaped relationship for the 5 dots - 1 dot and 3 dot - 1 dot, but those relationships are not significant.

library(tidyverse)

## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.2.1     ✓ purrr   0.3.3
## ✓ tibble  2.1.3     ✓ dplyr   0.8.3
## ✓ tidyr   1.0.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.4.0

## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(psych)

## 
## Attaching package: 'psych'

## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha

library(patchwork)
library(rmatio)
library(reshape2)

## 
## Attaching package: 'reshape2'

## The following object is masked from 'package:tidyr':
## 
##     smiths

load("data/behav.RData")
load("data/load_effects_DFR.RData")
load("data/split_groups_info.RData")

source("helper_fxns/load_EEG_data.R")
source("helper_fxns/corr_spectrum_to_indiv_diff.R")
source("helper_fxns/mutate_for_heatmap.R")
source("helper_fxns/split_ERPs_into_groups.R")
source("helper_fxns/create_TC_for_plot.R")
source("helper_fxns/split_spectrum.R")
source("helper_fxns/paired_freq_plot.R")
source("helper_fxns/select_period_average.R")
source("helper_fxns/split_into_groups.R")
source("helper_fxns/prep_split_for_bar_plots.R")
source("helper_fxns/average_electrodes.R")
source("helper_fxns/average_electrodes_spectrum.R")

temp <- read.mat("data/EEG/DFR/ERPS_cl15_midOccip_reformatted.mat")
ERPS_times_DFR <- temp[["all_times"]]
temp <- read.mat("data/EEG/DFR/ERSPS_cl15_midOccip_reformatted.mat")
ERSPS_times_DFR <- temp[["all_times"]]
ERSPS_freqs <- temp[["all_freqs"]]


ERPS_midOccip_DFR <- load_EEG_data("DFR", "ERPS_cl15_midOccip")
ERSPS_midOccip_DFR <- load_EEG_data("DFR", "ERSPS_cl15_midOccip")
ERSPS_Oz_DFR <- load_EEG_data("DFR", "ERSPS_Oz")

ERPS_Pz_DFR <- load_EEG_data("DFR", "ERPS_Pz")
ERPS_P1_DFR <- load_EEG_data("DFR", "ERPS_P1")
ERPS_P2_DFR <- load_EEG_data("DFR", "ERPS_P2")
ERPS_CP1_DFR <- load_EEG_data("DFR", "ERPS_CP1")
ERPS_CP2_DFR <- load_EEG_data("DFR", "ERPS_CP2")
ERPS_POz_DFR <- load_EEG_data("DFR", "ERPS_POz")

ERSPS_Fz_DFR <- load_EEG_data("DFR", "ERSPS_Fz")
ERSPS_F1_DFR <- load_EEG_data("DFR", "ERSPS_F1")
ERSPS_F2_DFR <- load_EEG_data("DFR", "ERSPS_F2")
ERSPS_FC1_DFR <- load_EEG_data("DFR", "ERSPS_FC1")
ERSPS_FC2_DFR <- load_EEG_data("DFR", "ERSPS_FC2")
ERSPS_AFz_DFR <- load_EEG_data("DFR", "ERSPS_AFz")


ERPS_O1_DFR <- load_EEG_data("DFR", "ERPS_O1")
ERPS_O2_DFR <- load_EEG_data("DFR", "ERPS_O2")
ERPS_PO7_DFR <- load_EEG_data("DFR", "ERPS_PO7")
ERPS_PO8_DFR <- load_EEG_data("DFR", "ERPS_PO8")


ERPS_Oavg_DFR <- average_electrodes(list(ERPS_O1_DFR, ERPS_O2_DFR))
ERPS_POavg_DFR <- average_electrodes(list(ERPS_PO7_DFR, ERPS_PO8_DFR))

ERPS_Pavg_DFR <- average_electrodes(list(ERPS_Pz_DFR,ERPS_P1_DFR, ERPS_P2_DFR,ERPS_POz_DFR,ERPS_CP1_DFR, ERPS_CP2_DFR))
ERSPS_Favg_DFR <- average_electrodes_spectrum(list(ERSPS_Fz_DFR, ERSPS_F1_DFR, ERSPS_F2_DFR, ERSPS_AFz_DFR, ERSPS_FC1_DFR, ERSPS_FC2_DFR))

#save(list=c("CDA","CDA_fMRI","ERPS_times_DFR", "ERPS_times_LCD", "ERSPS_times_DFR", "ERSPS_times_LCD", "ERSPS_freqs", "ERPS_midOccip_LCD", "ERSPS_midOccip_LCD", "ERSPS_Oz_LCD","ERSPS_O2_LCD","ERSPS_O1_LCD", "ERPS_Pz_LCD","ERPS_Fz_LCD", "ERPS_midOccip_DFR", "ERSPS_midOccip_DFR", "ERSPS_Oz_DFR","ERPS_Pz_DFR"), file="data/newEEG_data.RData")

rects <- data.frame(xstart=c(0,5500),xend=c(2500,7000),col = "gray")

ERP Data

avg_ERPs_for_plot <- list(
  midOccip_DFR = data.frame(high_load = apply(ERPS_midOccip_DFR[["high_load"]][,2:1922], 2, mean),
                            low_load = apply(ERPS_midOccip_DFR[["low_load"]][,2:1922], 2, mean),
                            load_effect = apply(ERPS_midOccip_DFR[["load_effect"]][,2:1922], 2, mean),
                            time =ERPS_times_DFR),
  Pavg_DFR = data.frame(high_load = apply(ERPS_Pavg_DFR[["high_load"]][,2:1922], 2, mean),
                        low_load = apply(ERPS_Pavg_DFR[["low_load"]][,2:1922], 2, mean),
                        load_effect = apply(ERPS_Pavg_DFR[["load_effect"]][,2:1922], 2, mean),
                        time =ERPS_times_DFR),
  Oavg_DFR = data.frame(high_load = apply(ERPS_Oavg_DFR[["high_load"]][,2:1922], 2, mean),
                        low_load = apply(ERPS_Oavg_DFR[["low_load"]][,2:1922], 2, mean),
                        load_effect = apply(ERPS_Oavg_DFR[["load_effect"]][,2:1922], 2, mean),
                        time =ERPS_times_DFR), 
  POavg_DFR = data.frame(high_load = apply(ERPS_POavg_DFR[["high_load"]][,2:1922], 2, mean),
                         low_load = apply(ERPS_POavg_DFR[["low_load"]][,2:1922], 2, mean),
                         load_effect = apply(ERPS_POavg_DFR[["load_effect"]][,2:1922], 2, mean),
                         time =ERPS_times_DFR)
  
)

plot_list <- list()

for (cluster in seq.int(1,length(avg_ERPs_for_plot))){
  plot_list[[names(avg_ERPs_for_plot)[cluster]]][["indiv_loads"]] <- ggplot(data = avg_ERPs_for_plot[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=high_load))+
    geom_line(aes(x=time,y=low_load),linetype="dotted")+
    ylab("ERP")+
    ggtitle("Low Load vs High Load")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
  plot_list[[names(avg_ERPs_for_plot)[cluster]]][["LE"]] <- ggplot(data = avg_ERPs_for_plot[[cluster]])+
    geom_line(aes(x=time,y=load_effect))+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    ylab("ERP")+
    ggtitle("Load Effect")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
}

avg_P <- (avg_ERPs_for_plot[["Pavg_DFR"]]$high_load + avg_ERPs_for_plot[["Pavg_DFR"]]$low_load)/2
avg_O <- (avg_ERPs_for_plot[["Oavg_DFR"]]$high_load + avg_ERPs_for_plot[["Oavg_DFR"]]$low_load)/2
avg_PO <- (avg_ERPs_for_plot[["POavg_DFR"]]$high_load + avg_ERPs_for_plot[["POavg_DFR"]]$low_load)/2

avg_plots <- data.frame(time = avg_ERPs_for_plot[["Pavg_DFR"]]$time, avg_P, avg_O, avg_PO)

ggplot(data = avg_plots)+
  geom_line(aes(x=time,y=avg_P))+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  ylab("ERP")+
  ggtitle("Average P electrodes ERP")+
  scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
  theme_classic()

ggplot(data = avg_plots)+
  geom_line(aes(x=time,y=avg_O))+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  ylab("ERP")+
  ggtitle("Average O electrodes ERP")+
  scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
  theme_classic()

ggplot(data = avg_plots)+
  geom_line(aes(x=time,y=avg_PO))+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  ylab("ERP")+
  ggtitle("Average PO electrodes ERP")+
  scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
  theme_classic()

cue_P3_range <- c(which(avg_plots$avg_P == max(avg_plots$avg_P[229:290]))-6 , which(avg_plots$avg_P == max(avg_plots$avg_P[229:290]))+6)


probe_P3_range <- c(which(avg_plots$avg_P == max(avg_plots$avg_P[1664:1680]))-6 , which(avg_plots$avg_P == max(avg_plots$avg_P[1664:1680]))+6)

cue_n170_range <- c(which(avg_plots$avg_O == min(avg_plots$avg_O[129:198]))-6, which(avg_plots$avg_O == min(avg_plots$avg_O[129:198]))+6)

probe_n170_range <- c(which(avg_plots$avg_O == min(avg_plots$avg_O[1550:1677]))-6,which(avg_plots$avg_O == min(avg_plots$avg_O[1550:1677]))+6)

All Subjects

During both the cue and probe, we see significant load effects in the N170 components in the average of the O1 and O2 electrodes and the PO7 and PO8 averages. We also see a load effect in the P3 in the parietal electrode average. Both of these show a larger amplitude in the low load than high load.

plot_list[["midOccip_DFR"]][["indiv_loads"]] + plot_list[["midOccip_DFR"]][["LE"]] +
  plot_annotation(title="midOccip - DFR")

plot_list[["Pavg_DFR"]][["indiv_loads"]] + plot_list[["Pavg_DFR"]][["LE"]] +
  plot_annotation(title="P avg - DFR")

plot_list[["Oavg_DFR"]][["indiv_loads"]] + plot_list[["Oavg_DFR"]][["LE"]] +
  plot_annotation(title="O avg - DFR")

plot_list[["POavg_DFR"]][["indiv_loads"]] + plot_list[["POavg_DFR"]][["LE"]] +
  plot_annotation(title="PO avg - DFR")

cue_average_midOccip_n170 <- select_period_average(ERPS_midOccip_DFR,avg_plots$time[cue_n170_range[1]],avg_plots$time[cue_n170_range[2]],ERPS_times_DFR)
cue_average_O_n170 <- select_period_average(ERPS_Oavg_DFR,avg_plots$time[cue_n170_range[1]],avg_plots$time[cue_n170_range[2]],ERPS_times_DFR)
cue_average_PO_n170 <- select_period_average(ERPS_POavg_DFR,avg_plots$time[cue_n170_range[1]],avg_plots$time[cue_n170_range[2]],ERPS_times_DFR)

cue_average_P3 <- select_period_average(ERPS_Pavg_DFR,avg_plots$time[cue_P3_range[1]],avg_plots$time[cue_P3_range[2]],ERPS_times_DFR)

probe_average_midOccip_n170 <- select_period_average(ERPS_midOccip_DFR,avg_plots$time[probe_n170_range[1]],avg_plots$time[probe_n170_range[2]],ERPS_times_DFR)
probe_average_O_n170 <- select_period_average(ERPS_Oavg_DFR,avg_plots$time[probe_n170_range[1]],avg_plots$time[probe_n170_range[2]],ERPS_times_DFR)
probe_average_PO_n170 <- select_period_average(ERPS_POavg_DFR,avg_plots$time[probe_n170_range[1]],avg_plots$time[probe_n170_range[2]],ERPS_times_DFR)

probe_average_P3 <- select_period_average(ERPS_Pavg_DFR,avg_plots$time[probe_P3_range[1]],avg_plots$time[probe_P3_range[2]],ERPS_times_DFR)

t.test(cue_average_midOccip_n170$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  cue_average_midOccip_n170$load_effect
## t = -1.0273, df = 177, p-value = 0.3057
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.08180836  0.02579485
## sample estimates:
##   mean of x 
## -0.02800676

t.test(cue_average_O_n170$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  cue_average_O_n170$load_effect
## t = -5.2605, df = 189, p-value = 3.871e-07
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.8283886 -0.3765553
## sample estimates:
## mean of x 
## -0.602472

t.test(cue_average_PO_n170$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  cue_average_PO_n170$load_effect
## t = -6.9534, df = 189, p-value = 5.642e-11
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.0704404 -0.5973194
## sample estimates:
##  mean of x 
## -0.8338799

t.test(cue_average_P3$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  cue_average_P3$load_effect
## t = -5.1679, df = 189, p-value = 5.988e-07
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.6709431 -0.3002420
## sample estimates:
##  mean of x 
## -0.4855926

t.test(probe_average_midOccip_n170$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  probe_average_midOccip_n170$load_effect
## t = 1.7801, df = 177, p-value = 0.07678
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.02215789  0.43009510
## sample estimates:
## mean of x 
## 0.2039686

t.test(probe_average_O_n170$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  probe_average_O_n170$load_effect
## t = 0.24254, df = 189, p-value = 0.8086
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.3504347  0.4486914
## sample estimates:
##  mean of x 
## 0.04912837

t.test(probe_average_PO_n170$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  probe_average_PO_n170$load_effect
## t = -0.28819, df = 189, p-value = 0.7735
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.8147653  0.6070440
## sample estimates:
##  mean of x 
## -0.1038607

t.test(probe_average_P3$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  probe_average_P3$load_effect
## t = -4.9474, df = 189, p-value = 1.657e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.7978863 -0.3430026
## sample estimates:
##  mean of x 
## -0.5704444

split_ERPS_midOccip_DFR <- split_ERPs_into_groups(ERPS_midOccip_DFR, WM_groups, ERPS_times_DFR) 
split_ERPS_Pavg_DFR <- split_ERPs_into_groups(ERPS_Pavg_DFR, WM_groups, ERPS_times_DFR) 
split_ERPS_Oavg_DFR <- split_ERPs_into_groups(ERPS_Oavg_DFR, WM_groups, ERPS_times_DFR) 
split_ERPS_POavg_DFR <- split_ERPs_into_groups(ERPS_POavg_DFR, WM_groups, ERPS_times_DFR) 


ERPs_for_plot <- list(mid_Occip_DFR = split_ERPS_midOccip_DFR,
                      Pavg_DFR = split_ERPS_Pavg_DFR, 
                      Oavg_DFR = split_ERPS_Oavg_DFR, 
                      POavg_DFR = split_ERPS_POavg_DFR)

split_ERPs_plot <- create_TC_for_plot(ERPs_for_plot)

Split by WM groups

ggplot(data = split_ERPs_plot[["mid_Occip_DFR"]][["long"]])+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  geom_line(data=split_ERPs_plot[["mid_Occip_DFR"]][["long"]] %>% filter(load=="high"),aes(x=Time,y=Mean,color=level)) +
  geom_line(data=split_ERPs_plot[["mid_Occip_DFR"]][["long"]] %>% filter(load=="low"),aes(x=Time,y=Mean,color=level),linetype="dotted")+
  scale_x_continuous(breaks = seq(-500,7000, by = 500))+
  ylab("Mean ERP")+
  ggtitle("mid_Occip_DFR")+
  theme_classic()

ggplot(data=split_ERPs_plot[["mid_Occip_DFR"]][["long"]])+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  geom_line(data=split_ERPs_plot[["mid_Occip_DFR"]][["long"]] %>% filter(load=="load_effect"),aes(x=Time,y=Mean,color=level)) +
  ylab("Mean Activity") +
  xlab("Time (ms)")+
  geom_ribbon(data=split_ERPs_plot[["mid_Occip_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="high"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="red")+
  geom_ribbon(data=split_ERPs_plot[["mid_Occip_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="med"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="green")+
  geom_ribbon(data=split_ERPs_plot[["mid_Occip_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="low"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="blue")+
  scale_x_continuous(breaks = seq(-500,7000, by = 500))+
  ylab("Mean ERP Load Effect")+
  ggtitle("mid_Occip_DFR")+
  theme_classic()

ggplot(data = split_ERPs_plot[["Pavg_DFR"]][["long"]])+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  geom_line(data=split_ERPs_plot[["Pavg_DFR"]][["long"]] %>% filter(load=="high"),aes(x=Time,y=Mean,color=level)) +
  geom_line(data=split_ERPs_plot[["Pavg_DFR"]][["long"]] %>% filter(load=="low"),aes(x=Time,y=Mean,color=level),linetype="dotted")+
  scale_x_continuous(breaks = seq(-500,7000, by = 500))+
  ylab("Mean ERP")+
  ggtitle("P avg DFR")+
  theme_classic()

ggplot(data=split_ERPs_plot[["Pavg_DFR"]][["long"]])+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  geom_line(data=split_ERPs_plot[["Pavg_DFR"]][["long"]] %>% filter(load=="load_effect"),aes(x=Time,y=Mean,color=level)) +
  ylab("Mean Activity") +
  xlab("Time (ms)")+
  geom_ribbon(data=split_ERPs_plot[["Pavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="high"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="red")+
  geom_ribbon(data=split_ERPs_plot[["Pavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="med"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="green")+
  geom_ribbon(data=split_ERPs_plot[["Pavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="low"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="blue")+
  scale_x_continuous(breaks = seq(-500,7000, by = 500))+
  ylab("Mean ERP Load Effect")+
  ggtitle("P avg DFR")+
  theme_classic()

ggplot(data = split_ERPs_plot[["Oavg_DFR"]][["long"]])+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  geom_line(data=split_ERPs_plot[["Oavg_DFR"]][["long"]] %>% filter(load=="high"),aes(x=Time,y=Mean,color=level)) +
  geom_line(data=split_ERPs_plot[["Oavg_DFR"]][["long"]] %>% filter(load=="low"),aes(x=Time,y=Mean,color=level),linetype="dotted")+
  scale_x_continuous(breaks = seq(-500,7000, by = 500))+
  ylab("Mean ERP")+
  ggtitle("O avg DFR")+
  theme_classic()

ggplot(data=split_ERPs_plot[["Oavg_DFR"]][["long"]])+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  geom_line(data=split_ERPs_plot[["Oavg_DFR"]][["long"]] %>% filter(load=="load_effect"),aes(x=Time,y=Mean,color=level)) +
  ylab("Mean Activity") +
  xlab("Time (ms)")+
  geom_ribbon(data=split_ERPs_plot[["Oavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="high"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="red")+
  geom_ribbon(data=split_ERPs_plot[["Oavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="med"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="green")+
  geom_ribbon(data=split_ERPs_plot[["Oavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="low"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="blue")+
  scale_x_continuous(breaks = seq(-500,7000, by = 500))+
  ylab("Mean ERP Load Effect")+
  ggtitle("O avg DFR")+
  theme_classic()

ggplot(data = split_ERPs_plot[["POavg_DFR"]][["long"]])+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  geom_line(data=split_ERPs_plot[["POavg_DFR"]][["long"]] %>% filter(load=="high"),aes(x=Time,y=Mean,color=level)) +
  geom_line(data=split_ERPs_plot[["POavg_DFR"]][["long"]] %>% filter(load=="low"),aes(x=Time,y=Mean,color=level),linetype="dotted")+
  scale_x_continuous(breaks = seq(-500,7000, by = 500))+
  ylab("Mean ERP")+
  ggtitle("PO avg DFR")+
  theme_classic()

ggplot(data=split_ERPs_plot[["POavg_DFR"]][["long"]])+
  geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
  geom_line(data=split_ERPs_plot[["POavg_DFR"]][["long"]] %>% filter(load=="load_effect"),aes(x=Time,y=Mean,color=level)) +
  ylab("Mean Activity") +
  xlab("Time (ms)")+
  geom_ribbon(data=split_ERPs_plot[["POavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="high"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="red")+
  geom_ribbon(data=split_ERPs_plot[["POavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="med"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="green")+
  geom_ribbon(data=split_ERPs_plot[["POavg_DFR"]][["long"]] %>% filter(load == "load_effect") %>% filter(level=="low"),aes(x=Time,ymin=SE_min, ymax=SE_max),alpha=.2,linetype=2,fill="blue")+
  scale_x_continuous(breaks = seq(-500,7000, by = 500))+
  ylab("Mean ERP Load Effect")+
  ggtitle("PO avg DFR")+
  theme_classic()

cue_P3_DFR_split <- split_into_groups(cue_average_P3, WM_groups)
cue_midOccip_n170_DFR_split <- split_into_groups(cue_average_midOccip_n170, WM_groups)
cue_Oavg_n170_DFR_split <- split_into_groups(cue_average_O_n170, WM_groups)
cue_POavg_n170_DFR_split <- split_into_groups(cue_average_PO_n170, WM_groups)


probe_P3_DFR_split <- split_into_groups(probe_average_P3, WM_groups)
probe_midOccip_n170_DFR_split <- split_into_groups(probe_average_midOccip_n170, WM_groups)
probe_Oavg_n170_DFR_split <- split_into_groups(probe_average_O_n170, WM_groups)
probe_POavg_n170_DFR_split <- split_into_groups(probe_average_PO_n170, WM_groups)


cue_P3_DFR_split_prepped <- prep_split_for_bar_plots(cue_P3_DFR_split)
cue_midOccip_n170_DFR_split_prepped <- prep_split_for_bar_plots(cue_midOccip_n170_DFR_split)
cue_Oavg_n170_DFR_split_prepped <- prep_split_for_bar_plots(cue_Oavg_n170_DFR_split)
cue_POavg_n170_DFR_split_prepped <- prep_split_for_bar_plots(cue_POavg_n170_DFR_split)


probe_P3_DFR_split_prepped <- prep_split_for_bar_plots(probe_P3_DFR_split)
probe_midOccip_n170_DFR_split_prepped <- prep_split_for_bar_plots(probe_midOccip_n170_DFR_split)
probe_Oavg_n170_DFR_split_prepped <- prep_split_for_bar_plots(probe_Oavg_n170_DFR_split)
probe_POavg_n170_DFR_split_prepped <- prep_split_for_bar_plots(probe_POavg_n170_DFR_split)

WM_to_merge <- WM_groups[["all"]][,c(1,12)]

cue_midOccip_n170_DFR_anova <- merge(cue_average_midOccip_n170,WM_to_merge, by="PTID")
probe_midOccip_n170_DFR_anova <- merge(probe_average_midOccip_n170,WM_to_merge, by="PTID")

cue_Oavg_n170_DFR_anova <- merge(cue_average_O_n170,WM_to_merge, by="PTID")
probe_Oavg_n170_DFR_anova <- merge(probe_average_O_n170,WM_to_merge, by="PTID")

cue_POavg_n170_DFR_anova <- merge(cue_average_PO_n170,WM_to_merge, by="PTID")
probe_POavg_n170_DFR_anova <- merge(probe_average_PO_n170,WM_to_merge, by="PTID")

cue_P3_DFR_anova <- merge(cue_average_P3,WM_to_merge, by="PTID")
probe_P3_DFR_anova <- merge(probe_average_P3,WM_to_merge, by="PTID")

High load

We see a linear increase in the P3 component at high load with capacity - there’s a larger amplitude with higher capacity.

We also do see the inverted U shaped relationship with the N170 component during cue for the occipital electrodes and component, but there’s high variance so it is not significant.

cue_P3_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> cue_P3_high

cue_midOccip_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> cue_midOccip_n170_high

cue_Oavg_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> cue_Oavg_n170_high

cue_POavg_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> cue_POavg_n170_high

probe_P3_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> probe_P3_high

probe_midOccip_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> probe_midOccip_n170_high

probe_Oavg_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> probe_Oavg_n170_high


probe_POavg_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> probe_POavg_n170_high

cue_P3_high + probe_P3_high+
  plot_annotation(title="P3 during high load during DFR")

cue_midOccip_n170_high + probe_midOccip_n170_high+
  plot_annotation(title="midOccip cluster n170 high load during DFR")

cue_Oavg_n170_high + probe_Oavg_n170_high+
  plot_annotation(title="Oavg n170 high load during DFR")

cue_POavg_n170_high + probe_POavg_n170_high+
  plot_annotation(title="POavg n170 high load during DFR")

cue_midOccip_n170_high.aov <- aov(high_load ~ level, data = cue_midOccip_n170_DFR_anova)
print(summary(cue_midOccip_n170_high.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.88  0.4385   0.598  0.551
## Residuals   148 108.48  0.7330

probe_midOccip_n170_high.aov <- aov(high_load ~ level, data = probe_midOccip_n170_DFR_anova)
print(summary(probe_midOccip_n170_high.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.86  0.4279   0.348  0.707
## Residuals   148 182.04  1.2300

cue_Oavg_n170_high.aov <- aov(high_load ~ level, data = cue_Oavg_n170_DFR_anova)
print(summary(cue_Oavg_n170_high.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   25.7   12.83   0.961  0.385
## Residuals   158 2110.5   13.36

probe_Oavg_n170_high.aov <- aov(high_load ~ level, data = probe_Oavg_n170_DFR_anova)
print(summary(probe_Oavg_n170_high.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## level         2   74.9   37.47   2.745 0.0673 .
## Residuals   158 2157.0   13.65                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

cue_POavg_n170_high.aov <- aov(high_load ~ level, data = cue_POavg_n170_DFR_anova)
print(summary(cue_POavg_n170_high.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   25.2   12.60   1.024  0.362
## Residuals   158 1944.7   12.31

probe_POavg_n170_high.aov <- aov(high_load ~ level, data = probe_POavg_n170_DFR_anova)
print(summary(probe_POavg_n170_high.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   44.3   22.17   1.609  0.203
## Residuals   158 2177.1   13.78

cue_P3_high.aov <- aov(high_load ~ level, data = cue_P3_DFR_anova)
print(summary(cue_P3_high.aov))

##              Df Sum Sq Mean Sq F value  Pr(>F)    
## level         2   73.2   36.58   9.394 0.00014 ***
## Residuals   158  615.2    3.89                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

probe_P3_high.aov <- aov(high_load ~ level, data = probe_P3_DFR_anova)
print(summary(probe_P3_high.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## level         2   44.4  22.182   4.461  0.013 *
## Residuals   158  785.7   4.973                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Load Effect

However, there is no difference in load effect

cue_P3_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> cue_P3_LE

cue_midOccip_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> cue_midOccip_n170_LE


cue_Oavg_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> cue_Oavg_n170_LE

cue_POavg_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> cue_POavg_n170_LE

probe_P3_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> probe_P3_LE

probe_midOccip_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> probe_midOccip_n170_LE

probe_Oavg_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> probe_Oavg_n170_LE

probe_POavg_n170_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> probe_POavg_n170_LE

cue_P3_LE + probe_P3_LE+
  plot_annotation(title="P3 Load Effect during DFR")

cue_midOccip_n170_LE + probe_midOccip_n170_LE+
  plot_annotation(title="midOccip cluster n170 Load Effect during DFR")

cue_Oavg_n170_LE + probe_Oavg_n170_LE+
  plot_annotation(title="Oavg n170 Load Effect during DFR")

cue_POavg_n170_LE + probe_POavg_n170_LE+
  plot_annotation(title="POavg n170 Load Effect during DFR")

cue_midOccip_n170_LE.aov <- aov(load_effect ~ level, data = cue_midOccip_n170_DFR_anova)
print(summary(cue_midOccip_n170_LE.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2  0.112  0.0560   0.439  0.646
## Residuals   148 18.882  0.1276

probe_midOccip_n170_LE.aov <- aov(load_effect ~ level, data = probe_midOccip_n170_DFR_anova)
print(summary(probe_midOccip_n170_LE.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    2.8   1.416   0.515  0.599
## Residuals   148  407.3   2.752

cue_Oavg_n170_LE.aov <- aov(load_effect ~ level, data = cue_Oavg_n170_DFR_anova)
print(summary(cue_Oavg_n170_LE.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    2.7   1.373   0.567  0.568
## Residuals   158  382.3   2.420

probe_Oavg_n170_LE.aov <- aov(load_effect ~ level, data = probe_Oavg_n170_DFR_anova)
print(summary(probe_Oavg_n170_LE.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    3.2   1.584   0.177  0.838
## Residuals   158 1412.4   8.939

cue_POavg_n170_LE.aov <- aov(load_effect ~ level, data = cue_POavg_n170_DFR_anova)
print(summary(cue_POavg_n170_LE.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    7.8   3.901   1.456  0.236
## Residuals   158  423.2   2.679

probe_POavg_n170_LE.aov <- aov(load_effect ~ level, data = probe_POavg_n170_DFR_anova)
print(summary(probe_POavg_n170_LE.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2     39   19.59   0.677   0.51
## Residuals   158   4571   28.93

cue_P3_LE.aov <- aov(load_effect ~ level, data = cue_P3_DFR_anova)
print(summary(cue_P3_LE.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.62  0.3077   0.181  0.835
## Residuals   158 268.74  1.7009

probe_P3_LE.aov <- aov(load_effect ~ level, data = probe_P3_DFR_anova)
print(summary(probe_P3_LE.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    0.3  0.1252   0.046  0.955
## Residuals   158  428.7  2.7134

ERSPs

temp_diff <- constructs_fMRI[,c(1,8)]

omnibus_span_corr <- list(
  ERSPS_midOccip_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_midOccip_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_midOccip_DFR[["PTID"]], times=ERSPS_times_DFR),
  ERSPS_Oz_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_Oz_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_Oz_DFR[["PTID"]], times=ERSPS_times_DFR), 
  
  ERSPS_Favg_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_Favg_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_Favg_DFR[["PTID"]], times=ERSPS_times_DFR) 
)

## Warning: `as_tibble.matrix()` requires a matrix with column names or a `.name_repair` argument. Using compatibility `.name_repair`.
## This warning is displayed once per session.

temp_diff <- p200_data[p200_data$PTID %in% constructs_fMRI$PTID, c(1,7)]

DFR_L3_acc_corr <- list(
  ERSPS_midOccip_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_midOccip_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_midOccip_DFR[["PTID"]], times=ERSPS_times_DFR),
  ERSPS_Oz_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_Oz_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_Oz_DFR[["PTID"]], times=ERSPS_times_DFR), 
  
  ERSPS_Favg_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_Favg_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_Favg_DFR[["PTID"]], times=ERSPS_times_DFR) 
)

temp_diff <- data.frame(PTID=constructs_fMRI$PTID, LE = p200_delay_DFR$DFR_Load3_Load1)

DFR_fMRI_LE_corr <- list(
  ERSPS_midOccip_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_midOccip_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_midOccip_DFR[["PTID"]], times=ERSPS_times_DFR),
  ERSPS_Oz_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_Oz_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_Oz_DFR[["PTID"]], times=ERSPS_times_DFR), 
  
  ERSPS_Favg_DFR = corr_spectrum_to_indiv_diff(indiv_diff = temp_diff, spectrum = ERSPS_Favg_DFR[["data"]][["load_effect"]], spec_PTID = ERSPS_Favg_DFR[["PTID"]], times=ERSPS_times_DFR) 
)

Average over all subs

avg_ERSPS_midOccip_DFR <- mutate_for_heatmap(apply(ERSPS_midOccip_DFR[["data"]][["load_effect"]],c(1,2),mean))
avg_ERSPS_Oz_DFR <- mutate_for_heatmap(apply(ERSPS_Oz_DFR[["data"]][["load_effect"]],c(1,2),mean))
avg_ERSPS_Favg_DFR <- mutate_for_heatmap(apply(ERSPS_Favg_DFR[["data"]][["load_effect"]],c(1,2),mean))


ggplot(data=avg_ERSPS_midOccip_DFR, aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("midOccip cluster during DFR")+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Power Load Effect")+
  theme_classic()

ggplot(data=avg_ERSPS_Oz_DFR, aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("Oz electrode during DFR")+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Power Load Effect")+
  theme_classic()

ggplot(data=avg_ERSPS_Favg_DFR, aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("F avg electrodes during DFR")+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Power Load Effect")+
  theme_classic()

Split by WMC

These plots give a good overview, but we will inspect each of the frequency bands more in depth later.

split_ERSPS_midOccip_DFR <- split_spectrum(ERSPS_midOccip_DFR, WM_groups)
split_ERSPS_Oz_DFR <- split_spectrum(ERSPS_Oz_DFR, WM_groups)
split_ERSPS_Favg_DFR <- split_spectrum(ERSPS_Favg_DFR, WM_groups)

split_ERSPS_midOccip_DFR_for_plot <- list()
split_ERSPS_Oz_DFR_for_plot <- list()
split_ERSPS_Favg_DFR_for_plot <- list()

split_ERSPS_midOccip_DFR_plots <- list()
split_ERSPS_Oz_DFR_plots <- list()
split_ERSPS_Favg_DFR_plots <- list()


for (group in c("high","med","low")){
  split_ERSPS_midOccip_DFR_for_plot[[group]] <- mutate_for_heatmap(split_ERSPS_midOccip_DFR[[group]][["load_effect"]])
  split_ERSPS_Oz_DFR_for_plot[[group]] <- mutate_for_heatmap(split_ERSPS_Oz_DFR[[group]][["load_effect"]])
  split_ERSPS_Favg_DFR_for_plot[[group]] <- mutate_for_heatmap(split_ERSPS_Favg_DFR[[group]][["load_effect"]])
  
  
  
  
  split_ERSPS_midOccip_DFR_plots[[group]] <- ggplot(data=split_ERSPS_midOccip_DFR_for_plot[[group]], aes(x=Y,y=X,fill=Z))+
    geom_tile()+
    ggtitle(paste0(group, " WMC in midOccip during DFR"))+
    geom_vline(linetype = "dotted",xintercept = 11)+
    geom_vline(linetype = "dotted", xintercept = 61)+
    geom_vline(linetype= "dotted", xintercept = 122)+
    scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
    scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
    scale_fill_gradient2(low="blue",mid="white",
                         high="red") +
    xlab("Time (ms)")+
    ylab("Frequency")+
    labs(fill="Power Load Effect")+
    theme_classic()
  
  split_ERSPS_Oz_DFR_plots[[group]] <- ggplot(data=split_ERSPS_Oz_DFR_for_plot[[group]], aes(x=Y,y=X,fill=Z))+
    geom_tile()+
    ggtitle(paste0(group, " WMC in Oz during DFR"))+
    geom_vline(linetype = "dotted",xintercept = 11)+
    geom_vline(linetype = "dotted", xintercept = 61)+
    geom_vline(linetype= "dotted", xintercept = 122)+
    scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
    scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
    scale_fill_gradient2(low="blue",mid="white",
                         high="red") +
    xlab("Time (ms)")+
    ylab("Frequency")+
    labs(fill="Power Load Effect")+
    theme_classic() 
  
  split_ERSPS_Favg_DFR_plots[[group]] <- ggplot(data=split_ERSPS_Favg_DFR_for_plot[[group]], aes(x=Y,y=X,fill=Z))+
    geom_tile()+
    ggtitle(paste0(group, " WMC in Favg during DFR"))+
    geom_vline(linetype = "dotted",xintercept = 11)+
    geom_vline(linetype = "dotted", xintercept = 61)+
    geom_vline(linetype= "dotted", xintercept = 122)+
    scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
    scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
    scale_fill_gradient2(low="blue",mid="white",
                         high="red") +
    xlab("Time (ms)")+
    ylab("Frequency")+
    labs(fill="Power Load Effect")+
    theme_classic() 
}

split_ERSPS_midOccip_DFR_plots[["high"]] 

split_ERSPS_midOccip_DFR_plots[["med"]] 

split_ERSPS_midOccip_DFR_plots[["low"]]

split_ERSPS_Oz_DFR_plots[["high"]]

split_ERSPS_Oz_DFR_plots[["med"]]

split_ERSPS_Oz_DFR_plots[["low"]]

split_ERSPS_Favg_DFR_plots[["high"]]

split_ERSPS_Favg_DFR_plots[["med"]]

split_ERSPS_Favg_DFR_plots[["low"]]

Correlation to Omnibus Span

Mid occipital cluster shows a positive correlation with span during the encoding and delay period in the beta and low gamma bands, and a negative correlation in the theta through alpha bands during the entire task, but most prominently right at the beginning of the probe period.

Oz shows a strong positive correlation with alpha power during encoding and what looks to be low gamma across the entire task. There also seems to be a negative correlation with beta across the entire task, with particularly strong negative correlations right at the end of the delay period/beginning of probe period.

F avg shows correlations during encoding with the theta band and low gamma/high beta that drift into the beginning of delay, and then a strong negative correlation right at the end of the delay period/beginning of the probe period.

ggplot(data=omnibus_span_corr[["ERSPS_midOccip_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("midOccip cluster during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=omnibus_span_corr[["ERSPS_Oz_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("Oz electrode during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=omnibus_span_corr[["ERSPS_Favg_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("F avg electrodes during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

Correlation to DFR high load accuracy

ggplot(data=DFR_L3_acc_corr[["ERSPS_midOccip_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("midOccip cluster during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=DFR_L3_acc_corr[["ERSPS_Oz_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("Oz electrode during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=DFR_L3_acc_corr[["ERSPS_Favg_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("F avg electrodes during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

Correlation to DFR fMRI load effect

There seems to be a strong correlation between fMRI load effect and power in the alpha and beta band across electrodes, particularly during encoding.

ggplot(data=DFR_fMRI_LE_corr[["ERSPS_midOccip_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("midOccip cluster during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=DFR_fMRI_LE_corr[["ERSPS_Oz_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("Oz electrode during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

ggplot(data=DFR_fMRI_LE_corr[["ERSPS_Favg_DFR"]][["plot"]],aes(x=Y,y=X,fill=Z))+
  geom_tile()+
  ggtitle("F avg electrodes during DFR")+
  geom_vline(linetype = "dotted",xintercept = 11)+
  geom_vline(linetype = "dotted", xintercept = 61)+
  geom_vline(linetype= "dotted", xintercept = 122)+
  scale_y_continuous(labels = function(x){return(round(ERSPS_freqs[x+1], digits=2))},breaks = seq(0, 76, by = 5))+
  scale_x_continuous(labels = function(x){return(ERSPS_times_DFR[x+1])}, breaks = seq(0,150, by = 25))+
  scale_fill_gradient2(low="blue",mid="white",
                       high="red") +
  xlab("Time (ms)")+
  ylab("Frequency")+
  labs(fill="Correlation")+
  theme_classic()

alpha <- list(
  
  ERSPS_midOccip_DFR = list(
    high_load = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["high_load"]][26:38,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["low_load"]][26:38,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["load_effect"]][26:38,,], c(2,3), mean)))),  
  
  ERSPS_Favg_DFR = list(
    high_load = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["high_load"]][26:38,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["low_load"]][26:38,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["load_effect"]][26:38,,], c(2,3), mean)))),  
  
  ERSPS_Oz_DFR = list(
    high_load = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["high_load"]][26:38,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["low_load"]][26:38,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["load_effect"]][26:38,,], c(2,3), mean))))
)

beta <- list(
  
  ERSPS_midOccip_DFR = list(
    high_load = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["high_load"]][37:45,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["low_load"]][37:45,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["load_effect"]][37:45,,], c(2,3), mean)))),  
  
  ERSPS_Favg_DFR = list(
    high_load = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["high_load"]][37:45,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["low_load"]][37:45,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["load_effect"]][37:45,,], c(2,3), mean)))),
  
  ERSPS_Oz_DFR = list(
    high_load = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["high_load"]][37:45,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["low_load"]][37:45,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["load_effect"]][37:45,,], c(2,3), mean))))
)

theta <- list(
  
  ERSPS_midOccip_DFR = list(
    high_load = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["high_load"]][1:25,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["low_load"]][1:25,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["load_effect"]][1:25,,], c(2,3), mean)))),  
  
  ERSPS_Favg_DFR = list(
    high_load = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["high_load"]][1:25,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["low_load"]][1:25,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["load_effect"]][1:25,,], c(2,3), mean)))),
  
  ERSPS_Oz_DFR = list(
    high_load = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["high_load"]][1:25,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["low_load"]][1:25,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["load_effect"]][1:25,,], c(2,3), mean))))
)

low_gamma <- list(
  
  ERSPS_midOccip_DFR = list(
    high_load = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["high_load"]][61:75,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["low_load"]][61:75,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_midOccip_DFR$PTID,t(apply(ERSPS_midOccip_DFR[["data"]][["load_effect"]][61:75,,], c(2,3), mean)))),  
  
  ERSPS_Favg_DFR = list(
    high_load = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["high_load"]][61:75,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["low_load"]][61:75,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_Favg_DFR$PTID,t(apply(ERSPS_Favg_DFR[["data"]][["load_effect"]][61:75,,], c(2,3), mean)))),
  
  ERSPS_Oz_DFR = list(
    high_load = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["high_load"]][61:75,,], c(2,3), mean))),
    low_load = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["low_load"]][61:75,,], c(2,3), mean))),
    load_effect = data.frame(PTID = ERSPS_Oz_DFR$PTID,t(apply(ERSPS_Oz_DFR[["data"]][["load_effect"]][61:75,,], c(2,3), mean))))
)

alpha_list <- list(
  ERSPS_midOccip_DFR = data.frame(
    high_load = apply(alpha[["ERSPS_midOccip_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(alpha[["ERSPS_midOccip_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(alpha[["ERSPS_midOccip_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR),  
  ERSPS_Oz_DFR = data.frame(
    high_load = apply(alpha[["ERSPS_Oz_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(alpha[["ERSPS_Oz_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(alpha[["ERSPS_Oz_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR),
  ERSPS_Favg_DFR = data.frame(
    high_load = apply(alpha[["ERSPS_Favg_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(alpha[["ERSPS_Favg_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(alpha[["ERSPS_Favg_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR)
)
beta_list <- list(
  ERSPS_midOccip_DFR = data.frame(
    high_load = apply(beta[["ERSPS_midOccip_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(beta[["ERSPS_midOccip_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(beta[["ERSPS_midOccip_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR),  
  ERSPS_Oz_DFR = data.frame(
    high_load = apply(beta[["ERSPS_Oz_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(beta[["ERSPS_Oz_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(beta[["ERSPS_Oz_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR),
  ERSPS_Favg_DFR = data.frame(
    high_load = apply(beta[["ERSPS_Favg_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(beta[["ERSPS_Favg_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(beta[["ERSPS_Favg_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR)
)


theta_list <- list(
  ERSPS_midOccip_DFR = data.frame(
    high_load = apply(theta[["ERSPS_midOccip_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(theta[["ERSPS_midOccip_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(theta[["ERSPS_midOccip_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR),  
  ERSPS_Oz_DFR = data.frame(
    high_load = apply(theta[["ERSPS_Oz_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(theta[["ERSPS_Oz_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(theta[["ERSPS_Oz_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR),
  ERSPS_Favg_DFR = data.frame(
    high_load = apply(theta[["ERSPS_Favg_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(theta[["ERSPS_Favg_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(theta[["ERSPS_Favg_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR)
)

low_gamma_list <- list(
  ERSPS_midOccip_DFR = data.frame(
    high_load = apply(low_gamma[["ERSPS_midOccip_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(low_gamma[["ERSPS_midOccip_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(low_gamma[["ERSPS_midOccip_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR),  
  ERSPS_Oz_DFR = data.frame(
    high_load = apply(low_gamma[["ERSPS_Oz_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(low_gamma[["ERSPS_Oz_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(low_gamma[["ERSPS_Oz_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR),
  ERSPS_Favg_DFR = data.frame(
    high_load = apply(low_gamma[["ERSPS_Favg_DFR"]][["high_load"]][,2:153], 2, mean),
    low_load = apply(low_gamma[["ERSPS_Favg_DFR"]][["low_load"]][,2:153], 2, mean),
    load_effect = apply(low_gamma[["ERSPS_Favg_DFR"]][["load_effect"]][,2:153], 2, mean),
    time = ERSPS_times_DFR)
)

alpha_plot_list <- list()
beta_plot_list <- list()
theta_plot_list <- list()
low_gamma_plot_list <- list()

for (cluster in seq.int(1,length(alpha))){
  alpha_plot_list[[names(alpha_list)[cluster]]][["indiv_loads"]] <- ggplot(data = alpha_list[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=high_load))+
    geom_line(aes(x=time,y=low_load),linetype="dotted")+
    ylab("alpha power")+
    ggtitle("Low Load vs High Load")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
  alpha_plot_list[[names(alpha_list)[cluster]]][["LE"]] <- ggplot(data = alpha_list[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=load_effect))+
    ylab("alpha power")+
    ggtitle("Load Effect")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
}

for (cluster in seq.int(1,length(beta))){
  beta_plot_list[[names(beta_list)[cluster]]][["indiv_loads"]] <- ggplot(data = beta_list[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=high_load))+
    geom_line(aes(x=time,y=low_load),linetype="dotted")+
    ylab("beta power")+
    ggtitle("Low Load vs High Load")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
  beta_plot_list[[names(beta_list)[cluster]]][["LE"]] <- ggplot(data = beta_list[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=load_effect))+
    ylab("beta power")+
    ggtitle("Load Effect")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
}

for (cluster in seq.int(1,length(theta))){
  theta_plot_list[[names(theta_list)[cluster]]][["indiv_loads"]] <- ggplot(data = theta_list[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=high_load))+
    geom_line(aes(x=time,y=low_load),linetype="dotted")+
    ylab("theta power")+
    ggtitle("Low Load vs High Load")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
  theta_plot_list[[names(theta_list)[cluster]]][["LE"]] <- ggplot(data = theta_list[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=load_effect))+
    ylab("theta power")+
    ggtitle("Load Effect")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
}


for (cluster in seq.int(1,length(low_gamma))){
  low_gamma_plot_list[[names(low_gamma_list)[cluster]]][["indiv_loads"]] <- ggplot(data = low_gamma_list[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=high_load))+
    geom_line(aes(x=time,y=low_load),linetype="dotted")+
    ylab("low_gamma power")+
    ggtitle("Low Load vs High Load")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
  low_gamma_plot_list[[names(low_gamma_list)[cluster]]][["LE"]] <- ggplot(data = low_gamma_list[[cluster]])+
    geom_rect(data=rects,aes(xmin=xstart, xmax=xend, ymin = -Inf, ymax=Inf,alpha =0.005),fill="grey",show.legend = FALSE)+
    geom_line(aes(x=time,y=load_effect))+
    ylab("low_gamma power")+
    ggtitle("Load Effect")+
    scale_x_continuous(breaks = seq(-500,7000, by = 1000))+
    theme_classic()
  
}

Alpha

All Subjects

When we look at the load effects over all subjects, we see significant load effects in the alpha band during the cue and probe period in Oz, mid Occipital cluster and F average electrodes and during the delay period in the Oz electrode. During the cue and probe periods, we see stronger power in the low load, but we see the opposite effect in the delay period.

alpha_plot_list[["ERSPS_midOccip_DFR"]][["indiv_loads"]] + alpha_plot_list[["ERSPS_midOccip_DFR"]][["LE"]] +
  plot_annotation(title="midOccip - alpha during DFR")

alpha_plot_list[["ERSPS_Oz_DFR"]][["indiv_loads"]] + alpha_plot_list[["ERSPS_Oz_DFR"]][["LE"]] +
  plot_annotation(title="Oz - alpha during DFR")

alpha_plot_list[["ERSPS_Favg_DFR"]][["indiv_loads"]] + alpha_plot_list[["ERSPS_Favg_DFR"]][["LE"]] +
  plot_annotation(title="F avg - alpha during DFR")

alpha_cue_average_midOccip <- select_period_average(alpha[["ERSPS_midOccip_DFR"]],0,2531.25,ERSPS_times_DFR)
alpha_cue_average_Oz <- select_period_average(alpha[["ERSPS_Oz_DFR"]],0,2531.25,ERSPS_times_DFR)
alpha_cue_average_Favg <- select_period_average(alpha[["ERSPS_Favg_DFR"]],0,2531.25,ERSPS_times_DFR)


alpha_delay_average_midOccip <- select_period_average(alpha[["ERSPS_midOccip_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)
alpha_delay_average_Oz <- select_period_average(alpha[["ERSPS_Oz_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)
alpha_delay_average_Favg <- select_period_average(alpha[["ERSPS_Favg_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)


alpha_probe_average_midOccip <- select_period_average(alpha[["ERSPS_midOccip_DFR"]],5511.71875,7000,ERSPS_times_DFR)
alpha_probe_average_Oz <- select_period_average(alpha[["ERSPS_Oz_DFR"]],5511.71875,7000,ERSPS_times_DFR)
alpha_probe_average_Favg <- select_period_average(alpha[["ERSPS_Favg_DFR"]],5511.71875,7000,ERSPS_times_DFR)

t.test(alpha_cue_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_cue_average_midOccip$load_effect
## t = -8.7699, df = 177, p-value = 1.454e-15
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.8545668 -0.5406145
## sample estimates:
##  mean of x 
## -0.6975907

t.test(alpha_cue_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_cue_average_Oz$load_effect
## t = -4.4787, df = 189, p-value = 1.298e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.4582411 -0.1780102
## sample estimates:
##  mean of x 
## -0.3181257

t.test(alpha_cue_average_Favg$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_cue_average_Favg$load_effect
## t = -4.4186, df = 189, p-value = 1.671e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -2.5142775 -0.9622461
## sample estimates:
## mean of x 
## -1.738262

t.test(alpha_delay_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_delay_average_midOccip$load_effect
## t = -0.1686, df = 177, p-value = 0.8663
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.1825836  0.1538409
## sample estimates:
##   mean of x 
## -0.01437135

t.test(alpha_delay_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_delay_average_Oz$load_effect
## t = 2.3001, df = 189, p-value = 0.02253
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.02393801 0.31225787
## sample estimates:
## mean of x 
## 0.1680979

t.test(alpha_delay_average_Favg$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_delay_average_Favg$load_effect
## t = -0.33285, df = 189, p-value = 0.7396
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.0635197  0.7564286
## sample estimates:
##  mean of x 
## -0.1535456

t.test(alpha_probe_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_probe_average_midOccip$load_effect
## t = -3.6404, df = 177, p-value = 0.0003572
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.6577849 -0.1953179
## sample estimates:
##  mean of x 
## -0.4265514

t.test(alpha_probe_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_probe_average_Oz$load_effect
## t = -5.941, df = 189, p-value = 1.338e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.7887671 -0.3955406
## sample estimates:
##  mean of x 
## -0.5921539

t.test(alpha_probe_average_Favg$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  alpha_probe_average_Favg$load_effect
## t = -5.2797, df = 189, p-value = 3.533e-07
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -6.133077 -2.796738
## sample estimates:
## mean of x 
## -4.464908

alphas_for_plot <- list()
betas_for_plot <- list()
thetas_for_plot <- list()
low_gammas_for_plot <- list()
load_list <- c("low_load", "high_load", "load_effect")
ERSPS_list <- list( 
  ERSPS_midOccip_DFR = ERSPS_midOccip_DFR, 
  ERSPS_Oz_DFR = ERSPS_Oz_DFR, 
  ERSPS_Favg_DFR = ERSPS_Favg_DFR)

for (cluster in seq.int(1,3)){
  temp <- list()
  for (load in seq.int(1,3)){
    temp[[load_list[load]]] <- alpha[[cluster]][[load]]
  }
  if (substr(names(alpha)[cluster],nchar(names(alpha)[cluster])-4,nchar(names(alpha)[cluster])-2) == "LCD"){
    alphas_for_plot[[names(ERSPS_list)[cluster]]] <- split_ERPs_into_groups(temp, WM_groups, ERSPS_times_LCD)
  } else{
    alphas_for_plot[[names(ERSPS_list)[cluster]]] <- split_ERPs_into_groups(temp, WM_groups, ERSPS_times_DFR)
  }
}

for (cluster in seq.int(1,3)){
  temp <- list()
  for (load in seq.int(1,3)){
    temp[[load_list[load]]] <- beta[[cluster]][[load]]
  }
  if (substr(names(beta)[cluster],nchar(names(beta)[cluster])-4,nchar(names(beta)[cluster])-2) == "LCD"){
    betas_for_plot[[names(ERSPS_list)[cluster]]] <- split_ERPs_into_groups(temp, WM_groups, ERSPS_times_LCD)
  } else{
    betas_for_plot[[names(ERSPS_list)[cluster]]] <- split_ERPs_into_groups(temp, WM_groups, ERSPS_times_DFR)
  }
}

for (cluster in seq.int(1,3)){
  temp <- list()
  for (load in seq.int(1,3)){
    temp[[load_list[load]]] <- theta[[cluster]][[load]]
  }
  if (substr(names(theta)[cluster],nchar(names(theta)[cluster])-4,nchar(names(theta)[cluster])-2) == "LCD"){
    thetas_for_plot[[names(ERSPS_list)[cluster]]] <- split_ERPs_into_groups(temp, WM_groups, ERSPS_times_LCD)
  } else{
    thetas_for_plot[[names(ERSPS_list)[cluster]]] <- split_ERPs_into_groups(temp, WM_groups, ERSPS_times_DFR)
  }
}

for (cluster in seq.int(1,3)){
  temp <- list()
  for (load in seq.int(1,3)){
    temp[[load_list[load]]] <- low_gamma[[cluster]][[load]]
  }
  if (substr(names(low_gamma)[cluster],nchar(names(low_gamma)[cluster])-4,nchar(names(low_gamma)[cluster])-2) == "LCD"){
    low_gammas_for_plot[[names(ERSPS_list)[cluster]]] <- split_ERPs_into_groups(temp, WM_groups, ERSPS_times_LCD)
  } else{
    low_gammas_for_plot[[names(ERSPS_list)[cluster]]] <- split_ERPs_into_groups(temp, WM_groups, ERSPS_times_DFR)
  }
}

By WM groups

When we split into the working memory capacity groups, we see differences across capacity in the mid occipital cluster during probe period, where we see a more negative load effect in the high capacity group vs the medium capacity group. We also see a significant difference between the low and high capacity subjects in the raw high load power, where the high capacity subjects show more negative power than the low load subjects.

We also can see an inverted U shape during the delay period in the Oz electrode, though it isn’t significant.

split_alphas_plot <- create_TC_for_plot(alphas_for_plot)
split_betas_plot <- create_TC_for_plot(betas_for_plot)
split_thetas_plot <- create_TC_for_plot(thetas_for_plot)
split_low_gammas_plot <- create_TC_for_plot(low_gammas_for_plot)

midOccip_DFR_alpha <- paired_freq_plot(split_alphas_plot[["ERSPS_midOccip_DFR"]][["long"]], c("high","low", "load_effect"))

Oz_DFR_alpha <- paired_freq_plot(split_alphas_plot[["ERSPS_Oz_DFR"]][["long"]], c("high","low", "load_effect"))
Favg_DFR_alpha <- paired_freq_plot(split_alphas_plot[["ERSPS_Favg_DFR"]][["long"]], c("high","low", "load_effect"))


midOccip_DFR_alpha[[1]] + midOccip_DFR_alpha[[2]]+
  plot_annotation(title="midOccip - alpha during DFR")+
  plot_layout(guides="collect")

Oz_DFR_alpha[[1]] + Oz_DFR_alpha[[2]]+
  plot_annotation(title="Oz - alpha during DFR")+
  plot_layout(guides="collect")

Favg_DFR_alpha[[1]] + Favg_DFR_alpha[[2]]+
  plot_annotation(title="F avg electrodes - alpha during DFR")+
  plot_layout(guides="collect")

alpha_cue_Oz_DFR_split <- split_into_groups(alpha_cue_average_Oz, WM_groups)
alpha_cue_midOccip_DFR_split <- split_into_groups(alpha_cue_average_midOccip, WM_groups)
alpha_cue_Favg_DFR_split <- split_into_groups(alpha_cue_average_Favg, WM_groups)

alpha_delay_Oz_DFR_split <- split_into_groups(alpha_delay_average_Oz, WM_groups)
alpha_delay_midOccip_DFR_split <- split_into_groups(alpha_delay_average_midOccip, WM_groups)
alpha_delay_Favg_DFR_split <- split_into_groups(alpha_delay_average_Favg, WM_groups)

alpha_probe_Oz_DFR_split <- split_into_groups(alpha_probe_average_Oz, WM_groups)
alpha_probe_midOccip_DFR_split <- split_into_groups(alpha_probe_average_midOccip, WM_groups)
alpha_probe_Favg_DFR_split <- split_into_groups(alpha_probe_average_Favg, WM_groups)


alpha_cue_Oz_DFR_split_prepped <- prep_split_for_bar_plots(alpha_cue_Oz_DFR_split)
alpha_cue_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(alpha_cue_midOccip_DFR_split)
alpha_cue_Favg_DFR_split_prepped <- prep_split_for_bar_plots(alpha_cue_Favg_DFR_split)


alpha_delay_Oz_DFR_split_prepped <- prep_split_for_bar_plots(alpha_delay_Oz_DFR_split)
alpha_delay_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(alpha_delay_midOccip_DFR_split)
alpha_delay_Favg_DFR_split_prepped <- prep_split_for_bar_plots(alpha_delay_Favg_DFR_split)


alpha_probe_Oz_DFR_split_prepped <- prep_split_for_bar_plots(alpha_probe_Oz_DFR_split)
alpha_probe_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(alpha_probe_midOccip_DFR_split)
alpha_probe_Favg_DFR_split_prepped <- prep_split_for_bar_plots(alpha_probe_Favg_DFR_split)

alpha_cue_midOccip_DFR_anova <- merge(alpha_cue_average_midOccip,WM_to_merge, by="PTID")
alpha_delay_midOccip_DFR_anova <- merge(alpha_delay_average_midOccip,WM_to_merge, by="PTID")
alpha_probe_midOccip_DFR_anova <- merge(alpha_probe_average_midOccip,WM_to_merge, by="PTID")

alpha_cue_Oz_DFR_anova <- merge(alpha_cue_average_Oz,WM_to_merge, by="PTID")
alpha_delay_Oz_DFR_anova <- merge(alpha_delay_average_Oz,WM_to_merge, by="PTID")
alpha_probe_Oz_DFR_anova <- merge(alpha_probe_average_Oz,WM_to_merge, by="PTID")

alpha_cue_Favg_DFR_anova <- merge(alpha_cue_average_Favg,WM_to_merge, by="PTID")
alpha_delay_Favg_DFR_anova <- merge(alpha_delay_average_Favg,WM_to_merge, by="PTID")
alpha_probe_Favg_DFR_anova <- merge(alpha_probe_average_Favg,WM_to_merge, by="PTID")

alpha_cue_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_cue_Oz_high

alpha_cue_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_cue_midOccip_high

alpha_cue_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_cue_Favg_high

alpha_delay_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_delay_Oz_high

alpha_delay_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_delay_midOccip_high

alpha_delay_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_delay_Favg_high

alpha_probe_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_probe_Oz_high

alpha_probe_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_probe_midOccip_high


alpha_probe_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_probe_Favg_high

alpha_cue_Oz_high + alpha_delay_Oz_high + alpha_probe_Oz_high+
  plot_annotation(title="Oz electrode alpha frequency in high load trials during DFR")

alpha_cue_Favg_high + alpha_delay_Favg_high + alpha_probe_Favg_high+
  plot_annotation(title="F avg electrodes alpha frequency high load trials during DFR")

alpha_cue_midOccip_high + alpha_delay_midOccip_high + alpha_probe_midOccip_high+
  plot_annotation(title="midOccip cluster alpha frequency high load trials during DFR")

alpha_cue_midOccip.aov <- aov(high_load ~ level, data = alpha_cue_midOccip_DFR_anova)
print("alpha - cue midOccip")

## [1] "alpha - cue midOccip"

print(summary(alpha_cue_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2      4   1.984   0.757  0.471
## Residuals   148    388   2.622

alpha_delay_midOccip.aov <- aov(high_load ~ level, data = alpha_delay_midOccip_DFR_anova)
print("alpha - delay midOccip")

## [1] "alpha - delay midOccip"

print(summary(alpha_delay_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    8.7   4.350   1.338  0.265
## Residuals   148  481.1   3.251

alpha_probe_midOccip.aov <- aov(high_load ~ level, data = alpha_probe_midOccip_DFR_anova)
print("alpha - probe midOccip")

## [1] "alpha - probe midOccip"

print(summary(alpha_probe_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## level         2   6.71   3.357   3.084 0.0487 *
## Residuals   148 161.11   1.089                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(TukeyHSD(alpha_probe_midOccip.aov))

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = high_load ~ level, data = alpha_probe_midOccip_DFR_anova)
## 
## $level
##                diff         lwr       upr     p adj
## low-high  0.5025709  0.01039164 0.9947501 0.0441720
## med-high  0.1361373 -0.35084905 0.6231237 0.7859313
## med-low  -0.3664335 -0.86558064 0.1327136 0.1945464

alpha_cue_Oz.aov <- aov(high_load ~ level, data = alpha_cue_Oz_DFR_anova)
print("alpha - cue Oz")

## [1] "alpha - cue Oz"

print(summary(alpha_cue_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   10.7   5.355   2.274  0.106
## Residuals   158  372.1   2.355

alpha_delay_Oz.aov <- aov(high_load ~ level, data = alpha_delay_Oz_DFR_anova)
print("alpha - delay Oz")

## [1] "alpha - delay Oz"

print(summary(alpha_delay_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    0.2  0.0926    0.03   0.97
## Residuals   158  484.9  3.0689

alpha_probe_Oz.aov <- aov(high_load ~ level, data = alpha_probe_Oz_DFR_anova)
print("alpha - probe Oz")

## [1] "alpha - probe Oz"

print(summary(alpha_probe_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.05  0.5256   0.321  0.726
## Residuals   158 258.37  1.6353

alpha_cue_Favg.aov <- aov(high_load ~ level, data = alpha_cue_Favg_DFR_anova)
print("alpha - cue Favg")

## [1] "alpha - cue Favg"

print(summary(alpha_cue_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   3.51   1.757   0.884  0.415
## Residuals   158 314.14   1.988

alpha_delay_Favg.aov <- aov(high_load ~ level, data = alpha_delay_Favg_DFR_anova)
print("alpha - delay Favg")

## [1] "alpha - delay Favg"

print(summary(alpha_delay_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2      1   0.520   0.152  0.859
## Residuals   158    540   3.418

alpha_probe_Favg.aov <- aov(high_load ~ level, data = alpha_probe_Favg_DFR_anova)
print("alpha - probe Favg")

## [1] "alpha - probe Favg"

print(summary(alpha_probe_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.01  0.5069   0.447   0.64
## Residuals   158 179.27  1.1346

alpha_cue_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_cue_Oz_LE

alpha_cue_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_cue_midOccip_LE

alpha_cue_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_cue_Favg_LE

alpha_delay_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_delay_Oz_LE

alpha_delay_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_delay_midOccip_LE

alpha_delay_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_delay_Favg_LE

alpha_probe_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_probe_Oz_LE

alpha_probe_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_probe_midOccip_LE


alpha_probe_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> alpha_probe_Favg_LE

alpha_cue_Oz_LE + alpha_delay_Oz_LE + alpha_probe_Oz_LE+
  plot_annotation(title="Oz electrode alpha frequency Load Effect during DFR")

alpha_cue_Favg_LE + alpha_delay_Favg_LE + alpha_probe_Favg_LE+
  plot_annotation(title="F avg electrodes alpha frequency Load Effect during DFR")

alpha_cue_midOccip_LE + alpha_delay_midOccip_LE + alpha_probe_midOccip_LE+
  plot_annotation(title="midOccip cluster alpha frequency Load Effect during DFR")

alpha_cue_midOccip.aov <- aov(load_effect ~ level, data = alpha_cue_midOccip_DFR_anova)
print("alpha - cue midOccip")

## [1] "alpha - cue midOccip"

print(summary(alpha_cue_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   2.99   1.494    1.31  0.273
## Residuals   148 168.85   1.141

alpha_delay_midOccip.aov <- aov(load_effect ~ level, data = alpha_delay_midOccip_DFR_anova)
print("alpha - delay midOccip")

## [1] "alpha - delay midOccip"

print(summary(alpha_delay_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   2.13   1.065   0.808  0.448
## Residuals   148 194.95   1.317

alpha_probe_midOccip.aov <- aov(load_effect ~ level, data = alpha_probe_midOccip_DFR_anova)
print("alpha - probe midOccip")

## [1] "alpha - probe midOccip"

print(summary(alpha_probe_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## level         2   16.0   8.019   3.094 0.0483 *
## Residuals   148  383.6   2.592                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(TukeyHSD(alpha_probe_midOccip.aov))

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = load_effect ~ level, data = alpha_probe_midOccip_DFR_anova)
## 
## $level
##               diff         lwr       upr     p adj
## low-high 0.6592005 -0.10025431 1.4186553 0.1029359
## med-high 0.7034747 -0.04796732 1.4549167 0.0717120
## med-low  0.0442742 -0.72593234 0.8144807 0.9898416

alpha_cue_Oz.aov <- aov(load_effect ~ level, data = alpha_cue_Oz_DFR_anova)
print("alpha - cue Oz")

## [1] "alpha - cue Oz"

print(summary(alpha_cue_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.65  0.8270   0.899  0.409
## Residuals   158 145.27  0.9194

alpha_delay_Oz.aov <- aov(load_effect ~ level, data = alpha_delay_Oz_DFR_anova)
print("alpha - delay Oz")

## [1] "alpha - delay Oz"

print(summary(alpha_delay_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   2.15   1.073   1.063  0.348
## Residuals   158 159.52   1.010

alpha_probe_Oz.aov <- aov(load_effect ~ level, data = alpha_probe_Oz_DFR_anova)
print("alpha - probe Oz")

## [1] "alpha - probe Oz"

print(summary(alpha_probe_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   7.76   3.881   2.218  0.112
## Residuals   158 276.44   1.750

alpha_cue_Favg.aov <- aov(load_effect ~ level, data = alpha_cue_Favg_DFR_anova)
print("alpha - cue Favg")

## [1] "alpha - cue Favg"

print(summary(alpha_cue_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2     23   11.51   0.404  0.669
## Residuals   158   4504   28.51

alpha_delay_Favg.aov <- aov(load_effect ~ level, data = alpha_delay_Favg_DFR_anova)
print("alpha - delay Favg")

## [1] "alpha - delay Favg"

print(summary(alpha_delay_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    107   53.28   1.258  0.287
## Residuals   158   6690   42.34

alpha_probe_Favg.aov <- aov(load_effect ~ level, data = alpha_probe_Favg_DFR_anova)
print("alpha - probe Favg")

## [1] "alpha - probe Favg"

print(summary(alpha_probe_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## level         2    704   351.9   2.402 0.0939 .
## Residuals   158  23149   146.5                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Beta

All Subjects

Overall, we see differences in load effects during cue and probe in all electrodes (low > high) and in delay during mid occipital and Oz (high > low).

beta_plot_list[["ERSPS_midOccip_DFR"]][["indiv_loads"]] + beta_plot_list[["ERSPS_midOccip_DFR"]][["LE"]] +
  plot_annotation(title="midOccip - beta during DFR")

beta_plot_list[["ERSPS_Oz_DFR"]][["indiv_loads"]] + beta_plot_list[["ERSPS_Oz_DFR"]][["LE"]] +
  plot_annotation(title="Oz - beta during DFR")

beta_plot_list[["ERSPS_Favg_DFR"]][["indiv_loads"]] + beta_plot_list[["ERSPS_Favg_DFR"]][["LE"]] +
  plot_annotation(title="F avg - beta during DFR")

beta_cue_average_midOccip <- select_period_average(beta[["ERSPS_midOccip_DFR"]],0,2531.25,ERSPS_times_DFR)
beta_cue_average_Oz <- select_period_average(beta[["ERSPS_Oz_DFR"]],0,2531.25,ERSPS_times_DFR)
beta_cue_average_Favg <- select_period_average(beta[["ERSPS_Favg_DFR"]],0,2531.25,ERSPS_times_DFR)


beta_delay_average_midOccip <- select_period_average(beta[["ERSPS_midOccip_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)
beta_delay_average_Oz <- select_period_average(beta[["ERSPS_Oz_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)
beta_delay_average_Favg <- select_period_average(beta[["ERSPS_Favg_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)


beta_probe_average_midOccip <- select_period_average(beta[["ERSPS_midOccip_DFR"]],5511.71875,7000,ERSPS_times_DFR)
beta_probe_average_Oz <- select_period_average(beta[["ERSPS_Oz_DFR"]],5511.71875,7000,ERSPS_times_DFR)
beta_probe_average_Favg <- select_period_average(beta[["ERSPS_Favg_DFR"]],5511.71875,7000,ERSPS_times_DFR)

t.test(beta_cue_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_cue_average_midOccip$load_effect
## t = -6.2428, df = 177, p-value = 3.092e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.5050572 -0.2624377
## sample estimates:
##  mean of x 
## -0.3837475

t.test(beta_cue_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_cue_average_Oz$load_effect
## t = -3.7009, df = 189, p-value = 0.0002817
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.3785452 -0.1153149
## sample estimates:
## mean of x 
##  -0.24693

t.test(beta_cue_average_Favg$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_cue_average_Favg$load_effect
## t = -2.6149, df = 189, p-value = 0.009646
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.6410224 -0.2297544
## sample estimates:
##  mean of x 
## -0.9353884

t.test(beta_delay_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_delay_average_midOccip$load_effect
## t = 2.1457, df = 177, p-value = 0.03326
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.01171767 0.28023834
## sample estimates:
## mean of x 
##  0.145978

t.test(beta_delay_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_delay_average_Oz$load_effect
## t = 6.7658, df = 189, p-value = 1.615e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.3191722 0.5818761
## sample estimates:
## mean of x 
## 0.4505242

t.test(beta_delay_average_Favg$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_delay_average_Favg$load_effect
## t = 1.3947, df = 189, p-value = 0.1647
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.2261368  1.3176975
## sample estimates:
## mean of x 
## 0.5457804

t.test(beta_probe_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_probe_average_midOccip$load_effect
## t = -5.0094, df = 177, p-value = 1.315e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.7269556 -0.3160592
## sample estimates:
##  mean of x 
## -0.5215074

t.test(beta_probe_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_probe_average_Oz$load_effect
## t = -6.1079, df = 189, p-value = 5.63e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.7840744 -0.4012606
## sample estimates:
##  mean of x 
## -0.5926675

t.test(beta_probe_average_Favg$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  beta_probe_average_Favg$load_effect
## t = -6.3629, df = 189, p-value = 1.459e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -6.796846 -3.579886
## sample estimates:
## mean of x 
## -5.188366

Split by WMC

There are, however, no differences across groups in beta for either the high load trials or load effects.

midOccip_DFR_beta <- paired_freq_plot(split_betas_plot[["ERSPS_midOccip_DFR"]][["long"]], c("high","low", "load_effect"))

Oz_DFR_beta <- paired_freq_plot(split_betas_plot[["ERSPS_Oz_DFR"]][["long"]], c("high","low", "load_effect"))
Favg_DFR_beta <- paired_freq_plot(split_betas_plot[["ERSPS_Favg_DFR"]][["long"]], c("high","low", "load_effect"))


midOccip_DFR_beta[[1]] + midOccip_DFR_beta[[2]]+
  plot_annotation(title="midOccip - beta during DFR")+
  plot_layout(guides="collect")

Oz_DFR_beta[[1]] + Oz_DFR_beta[[2]]+
  plot_annotation(title="Oz - beta during DFR")+
  plot_layout(guides="collect")

Favg_DFR_beta[[1]] + Favg_DFR_beta[[2]]+
  plot_annotation(title="F avg electrodes - beta during DFR")+
  plot_layout(guides="collect")

beta_cue_Oz_DFR_split <- split_into_groups(beta_cue_average_Oz, WM_groups)
beta_cue_midOccip_DFR_split <- split_into_groups(beta_cue_average_midOccip, WM_groups)
beta_cue_Favg_DFR_split <- split_into_groups(beta_cue_average_Favg, WM_groups)

beta_delay_Oz_DFR_split <- split_into_groups(beta_delay_average_Oz, WM_groups)
beta_delay_midOccip_DFR_split <- split_into_groups(beta_delay_average_midOccip, WM_groups)
beta_delay_Favg_DFR_split <- split_into_groups(beta_delay_average_Favg, WM_groups)

beta_probe_Oz_DFR_split <- split_into_groups(beta_probe_average_Oz, WM_groups)
beta_probe_midOccip_DFR_split <- split_into_groups(beta_probe_average_midOccip, WM_groups)
beta_probe_Favg_DFR_split <- split_into_groups(beta_probe_average_Favg, WM_groups)


beta_cue_Oz_DFR_split_prepped <- prep_split_for_bar_plots(beta_cue_Oz_DFR_split)
beta_cue_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(beta_cue_midOccip_DFR_split)
beta_cue_Favg_DFR_split_prepped <- prep_split_for_bar_plots(beta_cue_Favg_DFR_split)


beta_delay_Oz_DFR_split_prepped <- prep_split_for_bar_plots(beta_delay_Oz_DFR_split)
beta_delay_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(beta_delay_midOccip_DFR_split)
beta_delay_Favg_DFR_split_prepped <- prep_split_for_bar_plots(beta_delay_Favg_DFR_split)


beta_probe_Oz_DFR_split_prepped <- prep_split_for_bar_plots(beta_probe_Oz_DFR_split)
beta_probe_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(beta_probe_midOccip_DFR_split)
beta_probe_Favg_DFR_split_prepped <- prep_split_for_bar_plots(beta_probe_Favg_DFR_split)

beta_cue_midOccip_DFR_anova <- merge(beta_cue_average_midOccip,WM_to_merge, by="PTID")
beta_delay_midOccip_DFR_anova <- merge(beta_delay_average_midOccip,WM_to_merge, by="PTID")
beta_probe_midOccip_DFR_anova <- merge(beta_probe_average_midOccip,WM_to_merge, by="PTID")

beta_cue_Oz_DFR_anova <- merge(beta_cue_average_Oz,WM_to_merge, by="PTID")
beta_delay_Oz_DFR_anova <- merge(beta_delay_average_Oz,WM_to_merge, by="PTID")
beta_probe_Oz_DFR_anova <- merge(beta_probe_average_Oz,WM_to_merge, by="PTID")

beta_cue_Favg_DFR_anova <- merge(beta_cue_average_Favg,WM_to_merge, by="PTID")
beta_delay_Favg_DFR_anova <- merge(beta_delay_average_Favg,WM_to_merge, by="PTID")
beta_probe_Favg_DFR_anova <- merge(beta_probe_average_Favg,WM_to_merge, by="PTID")

beta_cue_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_cue_Oz_high

beta_cue_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_cue_midOccip_high

beta_cue_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_cue_Favg_high

beta_delay_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_delay_Oz_high

beta_delay_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_delay_midOccip_high

beta_delay_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_delay_Favg_high

beta_probe_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_probe_Oz_high

beta_probe_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_probe_midOccip_high


beta_probe_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_probe_Favg_high

beta_cue_Oz_high + beta_delay_Oz_high + beta_probe_Oz_high+
  plot_annotation(title="Oz electrode beta frequency in high load trials during DFR")

beta_cue_Favg_high + beta_delay_Favg_high + beta_probe_Favg_high+
  plot_annotation(title="F avg electrodes beta frequency high load trials during DFR")

beta_cue_midOccip_high + beta_delay_midOccip_high + beta_probe_midOccip_high+
  plot_annotation(title="midOccip cluster beta frequency high load trials during DFR")

beta_cue_midOccip.aov <- aov(high_load ~ level, data = beta_cue_midOccip_DFR_anova)
print("beta - cue midOccip")

## [1] "beta - cue midOccip"

print(summary(beta_cue_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    4.6   2.316   1.011  0.366
## Residuals   148  339.0   2.290

beta_delay_midOccip.aov <- aov(high_load ~ level, data = beta_delay_midOccip_DFR_anova)
print("beta - delay midOccip")

## [1] "beta - delay midOccip"

print(summary(beta_delay_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   4.06   2.029   1.294  0.277
## Residuals   148 232.03   1.568

beta_probe_midOccip.aov <- aov(high_load ~ level, data = beta_probe_midOccip_DFR_anova)
print("beta - probe midOccip")

## [1] "beta - probe midOccip"

print(summary(beta_probe_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   4.25   2.124   1.967  0.143
## Residuals   148 159.81   1.080

beta_cue_Oz.aov <- aov(high_load ~ level, data = beta_cue_Oz_DFR_anova)
print("beta - cue Oz")

## [1] "beta - cue Oz"

print(summary(beta_cue_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    9.1   4.552   2.184  0.116
## Residuals   158  329.4   2.085

beta_delay_Oz.aov <- aov(high_load ~ level, data = beta_delay_Oz_DFR_anova)
print("beta - delay Oz")

## [1] "beta - delay Oz"

print(summary(beta_delay_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.27  0.1365   0.075  0.928
## Residuals   158 287.05  1.8168

beta_probe_Oz.aov <- aov(high_load ~ level, data = beta_probe_Oz_DFR_anova)
print("beta - probe Oz")

## [1] "beta - probe Oz"

print(summary(beta_probe_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   2.34   1.171   0.635  0.531
## Residuals   158 291.37   1.844

beta_cue_Favg.aov <- aov(high_load ~ level, data = beta_cue_Favg_DFR_anova)
print("beta - cue Favg")

## [1] "beta - cue Favg"

print(summary(beta_cue_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    2.1   1.026   0.502  0.606
## Residuals   158  322.9   2.043

beta_delay_Favg.aov <- aov(high_load ~ level, data = beta_delay_Favg_DFR_anova)
print("beta - delay Favg")

## [1] "beta - delay Favg"

print(summary(beta_delay_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.16  0.0811   0.045  0.956
## Residuals   158 287.49  1.8195

beta_probe_Favg.aov <- aov(high_load ~ level, data = beta_probe_Favg_DFR_anova)
print("beta - probe Favg")

## [1] "beta - probe Favg"

print(summary(beta_probe_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   5.01   2.507   1.958  0.145
## Residuals   158 202.28   1.280

beta_cue_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_cue_Oz_LE

beta_cue_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_cue_midOccip_LE

beta_cue_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_cue_Favg_LE

beta_delay_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_delay_Oz_LE

beta_delay_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_delay_midOccip_LE

beta_delay_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_delay_Favg_LE

beta_probe_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_probe_Oz_LE

beta_probe_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_probe_midOccip_LE


beta_probe_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> beta_probe_Favg_LE

beta_cue_Oz_LE + beta_delay_Oz_LE + beta_probe_Oz_LE+
  plot_annotation(title="Oz electrode beta frequency Load Effect during DFR")

beta_cue_Favg_LE + beta_delay_Favg_LE + beta_probe_Favg_LE+
  plot_annotation(title="F avg electrodes beta frequency Load Effect during DFR")

beta_cue_midOccip_LE + beta_delay_midOccip_LE + beta_probe_midOccip_LE+
  plot_annotation(title="midOccip cluster beta frequency Load Effect during DFR")

beta_cue_midOccip.aov <- aov(load_effect ~ level, data = beta_cue_midOccip_DFR_anova)
print("beta - cue midOccip")

## [1] "beta - cue midOccip"

print(summary(beta_cue_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.15  0.0771   0.118  0.888
## Residuals   148  96.42  0.6515

beta_delay_midOccip.aov <- aov(load_effect ~ level, data = beta_delay_midOccip_DFR_anova)
print("beta - delay midOccip")

## [1] "beta - delay midOccip"

print(summary(beta_delay_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.09  0.0450   0.051   0.95
## Residuals   148 129.64  0.8759

beta_probe_midOccip.aov <- aov(load_effect ~ level, data = beta_probe_midOccip_DFR_anova)
print("beta - probe midOccip")

## [1] "beta - probe midOccip"

print(summary(beta_probe_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   8.04   4.020   1.886  0.155
## Residuals   148 315.42   2.131

beta_cue_Oz.aov <- aov(load_effect ~ level, data = beta_cue_Oz_DFR_anova)
print("beta - cue Oz")

## [1] "beta - cue Oz"

print(summary(beta_cue_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.46  0.2306   0.301   0.74
## Residuals   158 120.94  0.7655

beta_delay_Oz.aov <- aov(load_effect ~ level, data = beta_delay_Oz_DFR_anova)
print("beta - delay Oz")

## [1] "beta - delay Oz"

print(summary(beta_delay_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.59  0.7932   0.995  0.372
## Residuals   158 125.95  0.7971

beta_probe_Oz.aov <- aov(load_effect ~ level, data = beta_probe_Oz_DFR_anova)
print("beta - probe Oz")

## [1] "beta - probe Oz"

print(summary(beta_probe_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   2.04   1.020   0.589  0.556
## Residuals   158 273.47   1.731

beta_cue_Favg.aov <- aov(load_effect ~ level, data = beta_cue_Favg_DFR_anova)
print("beta - cue Favg")

## [1] "beta - cue Favg"

print(summary(beta_cue_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2     23   11.43   0.506  0.604
## Residuals   158   3571   22.60

beta_delay_Favg.aov <- aov(load_effect ~ level, data = beta_delay_Favg_DFR_anova)
print("beta - delay Favg")

## [1] "beta - delay Favg"

print(summary(beta_delay_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2     98   48.95   1.642  0.197
## Residuals   158   4711   29.81

beta_probe_Favg.aov <- aov(load_effect ~ level, data = beta_probe_Favg_DFR_anova)
print("beta - probe Favg")

## [1] "beta - probe Favg"

print(summary(beta_probe_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    499   249.4   1.786  0.171
## Residuals   158  22067   139.7

Theta

All Subjects

In the beta band, we can see significant load effects during the cue and delay periods in the Oz electrode (high load > low load), and in the probe period for the Oz electrode, mid Occipital cluster and F avg during the probe period (low > high).

theta_plot_list[["ERSPS_midOccip_DFR"]][["indiv_loads"]] + theta_plot_list[["ERSPS_midOccip_DFR"]][["LE"]] +
  plot_annotation(title="midOccip - theta during DFR")

theta_plot_list[["ERSPS_Oz_DFR"]][["indiv_loads"]] + theta_plot_list[["ERSPS_Oz_DFR"]][["LE"]] +
  plot_annotation(title="Oz - theta during DFR")

theta_plot_list[["ERSPS_Favg_DFR"]][["indiv_loads"]] + theta_plot_list[["ERSPS_Favg_DFR"]][["LE"]] +
  plot_annotation(title="F avg - theta during DFR")

theta_cue_average_midOccip <- select_period_average(theta[["ERSPS_midOccip_DFR"]],0,2531.25,ERSPS_times_DFR)
theta_cue_average_Oz <- select_period_average(theta[["ERSPS_Oz_DFR"]],0,2531.25,ERSPS_times_DFR)
theta_cue_average_Favg <- select_period_average(theta[["ERSPS_Favg_DFR"]],0,2531.25,ERSPS_times_DFR)


theta_delay_average_midOccip <- select_period_average(theta[["ERSPS_midOccip_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)
theta_delay_average_Oz <- select_period_average(theta[["ERSPS_Oz_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)
theta_delay_average_Favg <- select_period_average(theta[["ERSPS_Favg_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)


theta_probe_average_midOccip <- select_period_average(theta[["ERSPS_midOccip_DFR"]],5511.71875,7000,ERSPS_times_DFR)
theta_probe_average_Oz <- select_period_average(theta[["ERSPS_Oz_DFR"]],5511.71875,7000,ERSPS_times_DFR)
theta_probe_average_Favg <- select_period_average(theta[["ERSPS_Favg_DFR"]],5511.71875,7000,ERSPS_times_DFR)

t.test(theta_cue_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  theta_cue_average_midOccip$load_effect
## t = -0.4662, df = 177, p-value = 0.6416
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.1718366  0.1061636
## sample estimates:
##   mean of x 
## -0.03283653

t.test(theta_cue_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  theta_cue_average_Oz$load_effect
## t = 5.8367, df = 189, p-value = 2.28e-08
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.2575268 0.5204594
## sample estimates:
## mean of x 
## 0.3889931

t.test(theta_cue_average_Favg$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  theta_cue_average_Favg$load_effect
## t = -1.6775, df = 189, p-value = 0.09509
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.05204117  0.08504174
## sample estimates:
##  mean of x 
## -0.4834997

t.test(theta_delay_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  theta_delay_average_midOccip$load_effect
## t = -0.21113, df = 177, p-value = 0.833
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.1921162  0.1549821
## sample estimates:
##   mean of x 
## -0.01856705

t.test(theta_delay_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  theta_delay_average_Oz$load_effect
## t = 2.9431, df = 189, p-value = 0.003656
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.07886558 0.39945721
## sample estimates:
## mean of x 
## 0.2391614

t.test(theta_delay_average_Favg$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  theta_delay_average_Favg$load_effect
## t = -1.2367, df = 189, p-value = 0.2177
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.2826749  0.2941081
## sample estimates:
##  mean of x 
## -0.4942834

t.test(theta_probe_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  theta_probe_average_midOccip$load_effect
## t = -6.444, df = 177, p-value = 1.067e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.1390326 -0.6049424
## sample estimates:
##  mean of x 
## -0.8719875

t.test(theta_probe_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  theta_probe_average_Oz$load_effect
## t = -6.7285, df = 189, p-value = 1.986e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.0040782 -0.5488154
## sample estimates:
##  mean of x 
## -0.7764468

t.test(theta_probe_average_Favg$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  theta_probe_average_Favg$load_effect
## t = -6.7464, df = 189, p-value = 1.799e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -7.908032 -4.329772
## sample estimates:
## mean of x 
## -6.118902

Split by WMC

We see significant differences between WMC groups during the delay period in the mid Occipital cluster (negative linear relationship) and the cue period in the Oz electrode (larger load effects in the high capacity subjects).

In the raw power at high load, we also see a difference in the delay period in the mid occipital cluster, where there is significantly more negative power in the high capacity subjects than the low capacity subjects.

midOccip_DFR_theta <- paired_freq_plot(split_thetas_plot[["ERSPS_midOccip_DFR"]][["long"]], c("high","low", "load_effect"))

Oz_DFR_theta <- paired_freq_plot(split_thetas_plot[["ERSPS_Oz_DFR"]][["long"]], c("high","low", "load_effect"))
Favg_DFR_theta <- paired_freq_plot(split_thetas_plot[["ERSPS_Favg_DFR"]][["long"]], c("high","low", "load_effect"))


midOccip_DFR_theta[[1]] + midOccip_DFR_theta[[2]]+
  plot_annotation(title="midOccip - theta during DFR")+
  plot_layout(guides="collect")

Oz_DFR_theta[[1]] + Oz_DFR_theta[[2]]+
  plot_annotation(title="Oz - theta during DFR")+
  plot_layout(guides="collect")

Favg_DFR_theta[[1]] + Favg_DFR_theta[[2]]+
  plot_annotation(title="F avg electrodes - theta during DFR")+
  plot_layout(guides="collect")

theta_cue_Oz_DFR_split <- split_into_groups(theta_cue_average_Oz, WM_groups)
theta_cue_midOccip_DFR_split <- split_into_groups(theta_cue_average_midOccip, WM_groups)
theta_cue_Favg_DFR_split <- split_into_groups(theta_cue_average_Favg, WM_groups)

theta_delay_Oz_DFR_split <- split_into_groups(theta_delay_average_Oz, WM_groups)
theta_delay_midOccip_DFR_split <- split_into_groups(theta_delay_average_midOccip, WM_groups)
theta_delay_Favg_DFR_split <- split_into_groups(theta_delay_average_Favg, WM_groups)

theta_probe_Oz_DFR_split <- split_into_groups(theta_probe_average_Oz, WM_groups)
theta_probe_midOccip_DFR_split <- split_into_groups(theta_probe_average_midOccip, WM_groups)
theta_probe_Favg_DFR_split <- split_into_groups(theta_probe_average_Favg, WM_groups)


theta_cue_Oz_DFR_split_prepped <- prep_split_for_bar_plots(theta_cue_Oz_DFR_split)
theta_cue_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(theta_cue_midOccip_DFR_split)
theta_cue_Favg_DFR_split_prepped <- prep_split_for_bar_plots(theta_cue_Favg_DFR_split)


theta_delay_Oz_DFR_split_prepped <- prep_split_for_bar_plots(theta_delay_Oz_DFR_split)
theta_delay_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(theta_delay_midOccip_DFR_split)
theta_delay_Favg_DFR_split_prepped <- prep_split_for_bar_plots(theta_delay_Favg_DFR_split)


theta_probe_Oz_DFR_split_prepped <- prep_split_for_bar_plots(theta_probe_Oz_DFR_split)
theta_probe_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(theta_probe_midOccip_DFR_split)
theta_probe_Favg_DFR_split_prepped <- prep_split_for_bar_plots(theta_probe_Favg_DFR_split)

theta_cue_midOccip_DFR_anova <- merge(theta_cue_average_midOccip,WM_to_merge, by="PTID")
theta_delay_midOccip_DFR_anova <- merge(theta_delay_average_midOccip,WM_to_merge, by="PTID")
theta_probe_midOccip_DFR_anova <- merge(theta_probe_average_midOccip,WM_to_merge, by="PTID")

theta_cue_Oz_DFR_anova <- merge(theta_cue_average_Oz,WM_to_merge, by="PTID")
theta_delay_Oz_DFR_anova <- merge(theta_delay_average_Oz,WM_to_merge, by="PTID")
theta_probe_Oz_DFR_anova <- merge(theta_probe_average_Oz,WM_to_merge, by="PTID")

theta_cue_Favg_DFR_anova <- merge(theta_cue_average_Favg,WM_to_merge, by="PTID")
theta_delay_Favg_DFR_anova <- merge(theta_delay_average_Favg,WM_to_merge, by="PTID")
theta_probe_Favg_DFR_anova <- merge(theta_probe_average_Favg,WM_to_merge, by="PTID")

theta_cue_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_cue_Oz_high

theta_cue_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_cue_midOccip_high

theta_cue_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_cue_Favg_high

theta_delay_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_delay_Oz_high

theta_delay_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_delay_midOccip_high

theta_delay_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_delay_Favg_high

theta_probe_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_probe_Oz_high

theta_probe_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_probe_midOccip_high


theta_probe_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_probe_Favg_high

theta_cue_Oz_high + theta_delay_Oz_high + theta_probe_Oz_high+
  plot_annotation(title="Oz electrode theta frequency in high load trials during DFR")

theta_cue_Favg_high + theta_delay_Favg_high + theta_probe_Favg_high+
  plot_annotation(title="F avg electrodes theta frequency high load trials during DFR")

theta_cue_midOccip_high + theta_delay_midOccip_high + theta_probe_midOccip_high+
  plot_annotation(title="midOccip cluster theta frequency high load trials during DFR")

theta_cue_midOccip.aov <- aov(high_load ~ level, data = theta_cue_midOccip_DFR_anova)
print("theta - cue midOccip")

## [1] "theta - cue midOccip"

print(summary(theta_cue_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.92  0.4585   0.454  0.636
## Residuals   148 149.50  1.0101

theta_delay_midOccip.aov <- aov(high_load ~ level, data = theta_delay_midOccip_DFR_anova)
print("theta - delay midOccip")

## [1] "theta - delay midOccip"

print(summary(theta_delay_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## level         2  10.55   5.276   3.954 0.0212 *
## Residuals   148 197.47   1.334                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(TukeyHSD(theta_delay_midOccip.aov))

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = high_load ~ level, data = theta_delay_midOccip_DFR_anova)
## 
## $level
##                diff        lwr       upr     p adj
## low-high  0.6453087  0.1004122 1.1902052 0.0157139
## med-high  0.2634821 -0.2756653 0.8026296 0.4807732
## med-low  -0.3818266 -0.9344373 0.1707841 0.2338716

theta_probe_midOccip.aov <- aov(high_load ~ level, data = theta_probe_midOccip_DFR_anova)
print("theta - probe midOccip")

## [1] "theta - probe midOccip"

print(summary(theta_probe_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.77  0.8828   1.623  0.201
## Residuals   148  80.50  0.5439

theta_cue_Oz.aov <- aov(high_load ~ level, data = theta_cue_Oz_DFR_anova)
print("theta - cue Oz")

## [1] "theta - cue Oz"

print(summary(theta_cue_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   5.72   2.858   1.699  0.186
## Residuals   158 265.85   1.683

theta_delay_Oz.aov <- aov(high_load ~ level, data = theta_delay_Oz_DFR_anova)
print("theta - delay Oz")

## [1] "theta - delay Oz"

print(summary(theta_delay_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.94  0.9721   0.693  0.501
## Residuals   158 221.49  1.4018

theta_probe_Oz.aov <- aov(high_load ~ level, data = theta_probe_Oz_DFR_anova)
print("theta - probe Oz")

## [1] "theta - probe Oz"

print(summary(theta_probe_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   4.52   2.260   1.656  0.194
## Residuals   158 215.66   1.365

theta_cue_Favg.aov <- aov(high_load ~ level, data = theta_cue_Favg_DFR_anova)
print("theta - cue Favg")

## [1] "theta - cue Favg"

print(summary(theta_cue_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   2.42  1.2087   1.778  0.172
## Residuals   158 107.39  0.6797

theta_delay_Favg.aov <- aov(high_load ~ level, data = theta_delay_Favg_DFR_anova)
print("theta - delay Favg")

## [1] "theta - delay Favg"

print(summary(theta_delay_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    0.4  0.1992   0.227  0.797
## Residuals   158  138.8  0.8785

theta_probe_Favg.aov <- aov(high_load ~ level, data = theta_probe_Favg_DFR_anova)
print("theta - probe Favg")

## [1] "theta - probe Favg"

print(summary(theta_probe_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## level         2   4.27  2.1371   2.665 0.0727 .
## Residuals   158 126.69  0.8019                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

theta_cue_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_cue_Oz_LE

theta_cue_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_cue_midOccip_LE

theta_cue_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_cue_Favg_LE

theta_delay_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_delay_Oz_LE

theta_delay_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_delay_midOccip_LE

theta_delay_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_delay_Favg_LE

theta_probe_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_probe_Oz_LE

theta_probe_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_probe_midOccip_LE


theta_probe_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> theta_probe_Favg_LE

theta_cue_Oz_LE + theta_delay_Oz_LE + theta_probe_Oz_LE+
  plot_annotation(title="Oz electrode theta frequency Load Effect during DFR")

theta_cue_Favg_LE + theta_delay_Favg_LE + theta_probe_Favg_LE+
  plot_annotation(title="F avg electrodes theta frequency Load Effect during DFR")

theta_cue_midOccip_LE + theta_delay_midOccip_LE + theta_probe_midOccip_LE+
  plot_annotation(title="midOccip cluster theta frequency Load Effect during DFR")

theta_cue_midOccip.aov <- aov(load_effect ~ level, data = theta_cue_midOccip_DFR_anova)
print("theta - cue midOccip")

## [1] "theta - cue midOccip"

print(summary(theta_cue_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   2.43   1.215   1.377  0.255
## Residuals   148 130.54   0.882

theta_delay_midOccip.aov <- aov(load_effect ~ level, data = theta_delay_midOccip_DFR_anova)
print("theta - delay midOccip")

## [1] "theta - delay midOccip"

print(summary(theta_delay_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## level         2  10.25   5.123   3.735 0.0262 *
## Residuals   148 203.01   1.372                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

theta_probe_midOccip.aov <- aov(load_effect ~ level, data = theta_probe_midOccip_DFR_anova)
print("theta - probe midOccip")

## [1] "theta - probe midOccip"

print(summary(theta_probe_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   12.7   6.374   1.788  0.171
## Residuals   148  527.7   3.565

theta_cue_Oz.aov <- aov(load_effect ~ level, data = theta_cue_Oz_DFR_anova)
print("theta - cue Oz")

## [1] "theta - cue Oz"

print(summary(theta_cue_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)   
## level         2   7.62   3.811   4.965 0.0081 **
## Residuals   158 121.27   0.768                  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

theta_delay_Oz.aov <- aov(load_effect ~ level, data = theta_delay_Oz_DFR_anova)
print("theta - delay Oz")

## [1] "theta - delay Oz"

print(summary(theta_delay_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.87   0.434   0.371   0.69
## Residuals   158 184.61   1.168

theta_probe_Oz.aov <- aov(load_effect ~ level, data = theta_probe_Oz_DFR_anova)
print("theta - probe Oz")

## [1] "theta - probe Oz"

print(summary(theta_probe_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    0.4  0.1805   0.075  0.928
## Residuals   158  380.2  2.4064

theta_cue_Favg.aov <- aov(load_effect ~ level, data = theta_cue_Favg_DFR_anova)
print("theta - cue Favg")

## [1] "theta - cue Favg"

print(summary(theta_cue_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## level         2     84   42.02   2.587 0.0784 .
## Residuals   158   2567   16.24                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

theta_delay_Favg.aov <- aov(load_effect ~ level, data = theta_delay_Favg_DFR_anova)
print("theta - delay Favg")

## [1] "theta - delay Favg"

print(summary(theta_delay_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2     36   17.83   0.577  0.563
## Residuals   158   4882   30.90

theta_probe_Favg.aov <- aov(load_effect ~ level, data = theta_probe_Favg_DFR_anova)
print("theta - probe Favg")

## [1] "theta - probe Favg"

print(summary(theta_probe_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    114   56.81   0.314  0.731
## Residuals   158  28583  180.91

Low Gamma

All Subjects

Significant load effects during delay and probe in the mid Occipital cluster and Oz electrode (low > high) and F avg during probe (low > high).

low_gamma_plot_list[["ERSPS_midOccip_DFR"]][["indiv_loads"]] + low_gamma_plot_list[["ERSPS_midOccip_DFR"]][["LE"]] +
  plot_annotation(title="midOccip - low_gamma during DFR")

low_gamma_plot_list[["ERSPS_Oz_DFR"]][["indiv_loads"]] + low_gamma_plot_list[["ERSPS_Oz_DFR"]][["LE"]] +
  plot_annotation(title="Oz - low_gamma during DFR")

low_gamma_plot_list[["ERSPS_Favg_DFR"]][["indiv_loads"]] + low_gamma_plot_list[["ERSPS_Favg_DFR"]][["LE"]] +
  plot_annotation(title="F avg - low_gamma during DFR")

low_gamma_cue_average_midOccip <- select_period_average(low_gamma[["ERSPS_midOccip_DFR"]],0,2531.25,ERSPS_times_DFR)
low_gamma_cue_average_Oz <- select_period_average(low_gamma[["ERSPS_Oz_DFR"]],0,2531.25,ERSPS_times_DFR)
low_gamma_cue_average_Favg <- select_period_average(low_gamma[["ERSPS_Favg_DFR"]],0,2531.25,ERSPS_times_DFR)


low_gamma_delay_average_midOccip <- select_period_average(low_gamma[["ERSPS_midOccip_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)
low_gamma_delay_average_Oz <- select_period_average(low_gamma[["ERSPS_Oz_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)
low_gamma_delay_average_Favg <- select_period_average(low_gamma[["ERSPS_Favg_DFR"]],2531.25,5511.71875,ERSPS_times_DFR)


low_gamma_probe_average_midOccip <- select_period_average(low_gamma[["ERSPS_midOccip_DFR"]],5511.71875,7000,ERSPS_times_DFR)
low_gamma_probe_average_Oz <- select_period_average(low_gamma[["ERSPS_Oz_DFR"]],5511.71875,7000,ERSPS_times_DFR)
low_gamma_probe_average_Favg <- select_period_average(low_gamma[["ERSPS_Favg_DFR"]],5511.71875,7000,ERSPS_times_DFR)

t.test(low_gamma_cue_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_cue_average_midOccip$load_effect
## t = -0.14199, df = 177, p-value = 0.8872
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.06904871  0.05977946
## sample estimates:
##    mean of x 
## -0.004634626

t.test(low_gamma_cue_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_cue_average_Oz$load_effect
## t = -1.4566, df = 189, p-value = 0.1469
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.21352909  0.03212947
## sample estimates:
##   mean of x 
## -0.09069981

t.test(low_gamma_cue_average_Favg$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_cue_average_Favg$load_effect
## t = -0.70598, df = 189, p-value = 0.4811
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.8515969  0.4026952
## sample estimates:
##  mean of x 
## -0.2244509

t.test(low_gamma_delay_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_delay_average_midOccip$load_effect
## t = -2.4175, df = 177, p-value = 0.01664
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.19328508 -0.01954772
## sample estimates:
##  mean of x 
## -0.1064164

t.test(low_gamma_delay_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_delay_average_Oz$load_effect
## t = -2.0996, df = 189, p-value = 0.03709
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.295546088 -0.009216663
## sample estimates:
##  mean of x 
## -0.1523814

t.test(low_gamma_delay_average_Favg$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_delay_average_Favg$load_effect
## t = -1.5405, df = 189, p-value = 0.1251
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.6600909  0.2041998
## sample estimates:
##  mean of x 
## -0.7279456

t.test(low_gamma_probe_average_midOccip$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_probe_average_midOccip$load_effect
## t = -3.8094, df = 177, p-value = 0.0001919
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.6239510 -0.1980945
## sample estimates:
##  mean of x 
## -0.4110228

t.test(low_gamma_probe_average_Oz$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_probe_average_Oz$load_effect
## t = -3.7986, df = 189, p-value = 0.000196
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.6970924 -0.2205602
## sample estimates:
##  mean of x 
## -0.4588263

t.test(low_gamma_probe_average_Favg$load_effect,mu=0)

## 
##  One Sample t-test
## 
## data:  low_gamma_probe_average_Favg$load_effect
## t = -3.4732, df = 189, p-value = 0.0006376
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -4.731170 -1.303682
## sample estimates:
## mean of x 
## -3.017426

Split by WMC

In the low gamma band, the only load effect that shows any difference is Oz during cue, with a trend for an inverted U shape relationship.

In the raw power, however, we do see differences between groups in the cue and delay periods in the mid occipcital cluster and Oz electrodes. In the cue and delay for mid occipital and delay for Oz, we see a significant low > high difference, and in cue period for the mid occipital and cue and delay for Oz, we see a med > high differnece.

midOccip_DFR_low_gamma <- paired_freq_plot(split_low_gammas_plot[["ERSPS_midOccip_DFR"]][["long"]], c("high","low", "load_effect"))

Oz_DFR_low_gamma <- paired_freq_plot(split_low_gammas_plot[["ERSPS_Oz_DFR"]][["long"]], c("high","low", "load_effect"))
Favg_DFR_low_gamma <- paired_freq_plot(split_low_gammas_plot[["ERSPS_Favg_DFR"]][["long"]], c("high","low", "load_effect"))


midOccip_DFR_low_gamma[[1]] + midOccip_DFR_low_gamma[[2]]+
  plot_annotation(title="midOccip - low_gamma during DFR")+
  plot_layout(guides="collect")

Oz_DFR_low_gamma[[1]] + Oz_DFR_low_gamma[[2]]+
  plot_annotation(title="Oz - low_gamma during DFR")+
  plot_layout(guides="collect")

Favg_DFR_low_gamma[[1]] + Favg_DFR_low_gamma[[2]]+
  plot_annotation(title="F avg electrodes - low_gamma during DFR")+
  plot_layout(guides="collect")

low_gamma_cue_Oz_DFR_split <- split_into_groups(low_gamma_cue_average_Oz, WM_groups)
low_gamma_cue_midOccip_DFR_split <- split_into_groups(low_gamma_cue_average_midOccip, WM_groups)
low_gamma_cue_Favg_DFR_split <- split_into_groups(low_gamma_cue_average_Favg, WM_groups)

low_gamma_delay_Oz_DFR_split <- split_into_groups(low_gamma_delay_average_Oz, WM_groups)
low_gamma_delay_midOccip_DFR_split <- split_into_groups(low_gamma_delay_average_midOccip, WM_groups)
low_gamma_delay_Favg_DFR_split <- split_into_groups(low_gamma_delay_average_Favg, WM_groups)

low_gamma_probe_Oz_DFR_split <- split_into_groups(low_gamma_probe_average_Oz, WM_groups)
low_gamma_probe_midOccip_DFR_split <- split_into_groups(low_gamma_probe_average_midOccip, WM_groups)
low_gamma_probe_Favg_DFR_split <- split_into_groups(low_gamma_probe_average_Favg, WM_groups)


low_gamma_cue_Oz_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_cue_Oz_DFR_split)
low_gamma_cue_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_cue_midOccip_DFR_split)
low_gamma_cue_Favg_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_cue_Favg_DFR_split)


low_gamma_delay_Oz_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_delay_Oz_DFR_split)
low_gamma_delay_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_delay_midOccip_DFR_split)
low_gamma_delay_Favg_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_delay_Favg_DFR_split)


low_gamma_probe_Oz_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_probe_Oz_DFR_split)
low_gamma_probe_midOccip_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_probe_midOccip_DFR_split)
low_gamma_probe_Favg_DFR_split_prepped <- prep_split_for_bar_plots(low_gamma_probe_Favg_DFR_split)

low_gamma_cue_midOccip_DFR_anova <- merge(low_gamma_cue_average_midOccip,WM_to_merge, by="PTID")
low_gamma_delay_midOccip_DFR_anova <- merge(low_gamma_delay_average_midOccip,WM_to_merge, by="PTID")
low_gamma_probe_midOccip_DFR_anova <- merge(low_gamma_probe_average_midOccip,WM_to_merge, by="PTID")

low_gamma_cue_Oz_DFR_anova <- merge(low_gamma_cue_average_Oz,WM_to_merge, by="PTID")
low_gamma_delay_Oz_DFR_anova <- merge(low_gamma_delay_average_Oz,WM_to_merge, by="PTID")
low_gamma_probe_Oz_DFR_anova <- merge(low_gamma_probe_average_Oz,WM_to_merge, by="PTID")

low_gamma_cue_Favg_DFR_anova <- merge(low_gamma_cue_average_Favg,WM_to_merge, by="PTID")
low_gamma_delay_Favg_DFR_anova <- merge(low_gamma_delay_average_Favg,WM_to_merge, by="PTID")
low_gamma_probe_Favg_DFR_anova <- merge(low_gamma_probe_average_Favg,WM_to_merge, by="PTID")

low_gamma_cue_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_cue_Oz_high

low_gamma_cue_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_cue_midOccip_high

low_gamma_cue_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_cue_Favg_high

low_gamma_delay_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_delay_Oz_high

low_gamma_delay_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_delay_midOccip_high

low_gamma_delay_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_delay_Favg_high

low_gamma_probe_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_probe_Oz_high

low_gamma_probe_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_probe_midOccip_high


low_gamma_probe_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="high_load") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_probe_Favg_high

low_gamma_cue_Oz_high + low_gamma_delay_Oz_high + low_gamma_probe_Oz_high+
  plot_annotation(title="Oz electrode low_gamma frequency in high load trials during DFR")

low_gamma_cue_Favg_high + low_gamma_delay_Favg_high + low_gamma_probe_Favg_high+
  plot_annotation(title="F avg electrodes low_gamma frequency high load trials during DFR")

low_gamma_cue_midOccip_high + low_gamma_delay_midOccip_high + low_gamma_probe_midOccip_high+
  plot_annotation(title="midOccip cluster low_gamma frequency high load trials during DFR")

low_gamma_cue_midOccip.aov <- aov(high_load ~ level, data = low_gamma_cue_midOccip_DFR_anova)
print("low_gamma - cue midOccip")

## [1] "low_gamma - cue midOccip"

print(summary(low_gamma_cue_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## level         2  1.833  0.9167    4.34 0.0147 *
## Residuals   148 31.257  0.2112                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(TukeyHSD(low_gamma_cue_midOccip.aov))

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = high_load ~ level, data = low_gamma_cue_midOccip_DFR_anova)
## 
## $level
##                 diff         lwr       upr     p adj
## low-high  0.23611702  0.01932486 0.4529092 0.0291757
## med-high  0.22549227  0.01098742 0.4399971 0.0368705
## med-low  -0.01062475 -0.23048608 0.2092366 0.9928098

low_gamma_delay_midOccip.aov <- aov(high_load ~ level, data = low_gamma_delay_midOccip_DFR_anova)
print("low_gamma - delay midOccip")

## [1] "low_gamma - delay midOccip"

print(summary(low_gamma_delay_midOccip.aov))

##              Df Sum Sq Mean Sq F value  Pr(>F)   
## level         2  1.377  0.6885   5.005 0.00788 **
## Residuals   148 20.358  0.1376                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(TukeyHSD(low_gamma_delay_midOccip.aov))

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = high_load ~ level, data = low_gamma_delay_midOccip_DFR_anova)
## 
## $level
##                diff         lwr        upr     p adj
## low-high  0.2332313  0.05827371 0.40818888 0.0054711
## med-high  0.1251381 -0.04797354 0.29824977 0.2043084
## med-low  -0.1080932 -0.28552768 0.06934131 0.3219788

low_gamma_probe_midOccip.aov <- aov(high_load ~ level, data = low_gamma_probe_midOccip_DFR_anova)
print("low_gamma - probe midOccip")

## [1] "low_gamma - probe midOccip"

print(summary(low_gamma_probe_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2  0.735  0.3675   2.052  0.132
## Residuals   148 26.500  0.1791

low_gamma_cue_Oz.aov <- aov(high_load ~ level, data = low_gamma_cue_Oz_DFR_anova)
print("low_gamma - cue Oz")

## [1] "low_gamma - cue Oz"

print(summary(low_gamma_cue_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## level         2  12.64   6.318   3.994 0.0203 *
## Residuals   158 249.95   1.582                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(TukeyHSD(low_gamma_cue_Oz.aov))

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = high_load ~ level, data = low_gamma_cue_Oz_DFR_anova)
## 
## $level
##               diff         lwr       upr     p adj
## low-high 0.4878716 -0.08772857 1.0634718 0.1142282
## med-high 0.6552395  0.08514466 1.2253343 0.0197941
## med-low  0.1673678 -0.41081666 0.7455523 0.7726439

low_gamma_delay_Oz.aov <- aov(high_load ~ level, data = low_gamma_delay_Oz_DFR_anova)
print("low_gamma - delay Oz")

## [1] "low_gamma - delay Oz"

print(summary(low_gamma_delay_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## level         2   9.67   4.834   4.558 0.0119 *
## Residuals   158 167.55   1.060                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

print(TukeyHSD(low_gamma_delay_Oz.aov))

##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = high_load ~ level, data = low_gamma_delay_Oz_DFR_anova)
## 
## $level
##                 diff         lwr       upr     p adj
## low-high  0.53339700  0.06213393 1.0046601 0.0221939
## med-high  0.49896157  0.03220594 0.9657172 0.0330488
## med-low  -0.03443543 -0.50781435 0.4389435 0.9838034

low_gamma_probe_Oz.aov <- aov(high_load ~ level, data = low_gamma_probe_Oz_DFR_anova)
print("low_gamma - probe Oz")

## [1] "low_gamma - probe Oz"

print(summary(low_gamma_probe_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.88  0.4396    0.29  0.748
## Residuals   158 239.22  1.5140

low_gamma_cue_Favg.aov <- aov(high_load ~ level, data = low_gamma_cue_Favg_DFR_anova)
print("low_gamma - cue Favg")

## [1] "low_gamma - cue Favg"

print(summary(low_gamma_cue_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   3.63   1.816   1.517  0.223
## Residuals   158 189.17   1.197

low_gamma_delay_Favg.aov <- aov(high_load ~ level, data = low_gamma_delay_Favg_DFR_anova)
print("low_gamma - delay Favg")

## [1] "low_gamma - delay Favg"

print(summary(low_gamma_delay_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   3.84   1.921    1.82  0.165
## Residuals   158 166.76   1.056

low_gamma_probe_Favg.aov <- aov(high_load ~ level, data = low_gamma_probe_Favg_DFR_anova)
print("low_gamma - probe Favg")

## [1] "low_gamma - probe Favg"

print(summary(low_gamma_probe_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.19  0.5934     0.5  0.608
## Residuals   158 187.66  1.1877

low_gamma_cue_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_cue_Oz_LE

low_gamma_cue_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_cue_midOccip_LE

low_gamma_cue_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Cue")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_cue_Favg_LE

low_gamma_delay_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_delay_Oz_LE

low_gamma_delay_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_delay_midOccip_LE

low_gamma_delay_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Delay")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_delay_Favg_LE

low_gamma_probe_Oz_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_probe_Oz_LE

low_gamma_probe_midOccip_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_probe_midOccip_LE


low_gamma_probe_Favg_DFR_split_prepped[["melt_data"]] %>% 
  filter(variable=="load_effect") %>% 
  mutate(level = factor(level, levels =c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
  ggtitle("Probe")+ 
  xlab("WMC")+
  ylab("ERP amplitude")+
  theme_classic()+
  theme(aspect.ratio = 1) -> low_gamma_probe_Favg_LE

low_gamma_cue_Oz_LE + low_gamma_delay_Oz_LE + low_gamma_probe_Oz_LE+
  plot_annotation(title="Oz electrode low_gamma frequency Load Effect during DFR")

low_gamma_cue_Favg_LE + low_gamma_delay_Favg_LE + low_gamma_probe_Favg_LE+
  plot_annotation(title="F avg electrodes low_gamma frequency Load Effect during DFR")

low_gamma_cue_midOccip_LE + low_gamma_delay_midOccip_LE + low_gamma_probe_midOccip_LE+
  plot_annotation(title="midOccip cluster low_gamma frequency Load Effect during DFR")

low_gamma_cue_midOccip.aov <- aov(load_effect ~ level, data = low_gamma_cue_midOccip_DFR_anova)
print("low_gamma - cue midOccip")

## [1] "low_gamma - cue midOccip"

print(summary(low_gamma_cue_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2  0.075 0.03768   0.222  0.801
## Residuals   148 25.132 0.16981

low_gamma_delay_midOccip.aov <- aov(load_effect ~ level, data = low_gamma_delay_midOccip_DFR_anova)
print("low_gamma - delay midOccip")

## [1] "low_gamma - delay midOccip"

print(summary(low_gamma_delay_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.37  0.1856   0.577  0.563
## Residuals   148  47.57  0.3214

low_gamma_probe_midOccip.aov <- aov(load_effect ~ level, data = low_gamma_probe_midOccip_DFR_anova)
print("low_gamma - probe midOccip")

## [1] "low_gamma - probe midOccip"

print(summary(low_gamma_probe_midOccip.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    4.5   2.255   0.964  0.384
## Residuals   148  346.3   2.340

low_gamma_cue_Oz.aov <- aov(load_effect ~ level, data = low_gamma_cue_Oz_DFR_anova)
print("low_gamma - cue Oz")

## [1] "low_gamma - cue Oz"

print(summary(low_gamma_cue_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)  
## level         2   3.58  1.7882   2.724 0.0687 .
## Residuals   158 103.73  0.6565                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

low_gamma_delay_Oz.aov <- aov(load_effect ~ level, data = low_gamma_delay_Oz_DFR_anova)
print("low_gamma - delay Oz")

## [1] "low_gamma - delay Oz"

print(summary(low_gamma_delay_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.64  0.3200   0.386   0.68
## Residuals   158 130.95  0.8288

low_gamma_probe_Oz.aov <- aov(load_effect ~ level, data = low_gamma_probe_Oz_DFR_anova)
print("low_gamma - probe Oz")

## [1] "low_gamma - probe Oz"

print(summary(low_gamma_probe_Oz.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    1.9  0.9627   0.335  0.716
## Residuals   158  454.3  2.8753

low_gamma_cue_Favg.aov <- aov(load_effect ~ level, data = low_gamma_cue_Favg_DFR_anova)
print("low_gamma - cue Favg")

## [1] "low_gamma - cue Favg"

print(summary(low_gamma_cue_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   63.1   31.57   2.026  0.135
## Residuals   158 2462.7   15.59

low_gamma_delay_Favg.aov <- aov(load_effect ~ level, data = low_gamma_delay_Favg_DFR_anova)
print("low_gamma - delay Favg")

## [1] "low_gamma - delay Favg"

print(summary(low_gamma_delay_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2     81   40.74   1.167  0.314
## Residuals   158   5515   34.90

low_gamma_probe_Favg.aov <- aov(load_effect ~ level, data = low_gamma_probe_Favg_DFR_anova)
print("low_gamma - probe Favg")

## [1] "low_gamma - probe Favg"

print(summary(low_gamma_probe_Favg.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    250   124.8   0.811  0.446
## Residuals   158  24298   153.8

PAC

The only significant PAC measure across WMC groups is the theta - beta coupling in the Fz electrode, with an inverted U shaped relationship.

PAC_midOccip <- list(high_load=read_delim("data/EEG/DFR/PAC/PAC_cl15_midOccip_maint_design1_highload.txt", delim=","), low_load = read_delim("data/EEG/DFR/PAC/PAC_cl15_midOccip_maint_design1_lowload.txt", delim=",") ) 

PAC_midFront <- list(high_load=read_delim("data/EEG/DFR/PAC/PAC_cl5_midFront_maint_design1_highload.txt", delim=","), low_load = read_delim("data/EEG/DFR/PAC/PAC_cl5_midFront_maint_design1_lowload.txt", delim=",") ) 

PAC_Oz <- list(high_load=read_delim("data/EEG/DFR/PAC/PAC_Oz_maint_design1_highload.txt", delim=","), low_load = read_delim("data/EEG/DFR/PAC/PAC_Oz_maint_design1_lowload.txt", delim=",") ) 

PAC_O1 <- list(high_load=read_delim("data/EEG/DFR/PAC/PAC_O1_maint_design1_highload.txt", delim=","), low_load = read_delim("data/EEG/DFR/PAC/PAC_O1_maint_design1_lowload.txt", delim=",") ) 

PAC_O2 <- list(high_load=read_delim("data/EEG/DFR/PAC/PAC_O2_maint_design1_highload.txt", delim=","), low_load = read_delim("data/EEG/DFR/PAC/PAC_O2_maint_design1_lowload.txt", delim=",") ) 

PAC_Fz <- list(high_load=read_delim("data/EEG/DFR/PAC/PAC_Fz_maint_design1_highload.txt", delim=","), low_load = read_delim("data/EEG/DFR/PAC/PAC_Fz_maint_design1_lowload.txt", delim=",") ) 

PAC_list <- list(midOccip = PAC_midOccip, midFront = PAC_midFront, Oz = PAC_Oz, O1 = PAC_O1, O2 = PAC_O2, Fz = PAC_Fz)

for (PAC in seq.int(1,length(PAC_list))){
  for (level in seq.int(1,2)){
    colnames(PAC_list[[PAC]][[level]]) <- c("PTID", "thetP_hbetA_mnt", "thetP_lgamA_mnt", "thetP_hgamA_mnt", "thetP_ugamA_mnt", "alphP_hbetA_mnt", "alphP_lgamA_mnt", "alphP_hgamA_mnt", "alphP_ugamA_mnt")
    PAC_list[[PAC]][[level]] <- PAC_list[[PAC]][[level]][,1:9]
  }
  
  PAC_list[[PAC]][["load_effect"]] <- PAC_list[[PAC]][["high_load"]] - PAC_list[[PAC]][["low_load"]]
  PAC_list[[PAC]][["load_effect"]]$PTID <- PAC_list[[PAC]][["high_load"]]$PTID
  
  PAC_list[[PAC]][["load_effect"]] <- merge(PAC_list[[PAC]][["load_effect"]], constructs_fMRI, by = "PTID")
  PAC_list[[PAC]][["load_effect"]] <- merge(PAC_list[[PAC]][["load_effect"]], WM_to_merge, by = "PTID")
  
  PAC_list[[PAC]][["high_load"]] <- merge(PAC_list[[PAC]][["high_load"]], constructs_fMRI, by = "PTID")
  PAC_list[[PAC]][["high_load"]] <- merge(PAC_list[[PAC]][["high_load"]], WM_to_merge, by = "PTID")
  
}

Theta

High Load

for (PAC in seq.int(1,length(PAC_list))){
  
  print(ggplot(data = PAC_list[[PAC]][["high_load"]], aes(x = omnibus_span_no_DFR, y = thetP_hbetA_mnt))+
          geom_point()+
          stat_smooth(method="lm")+
          ggtitle(names(PAC_list)[PAC]))
  
  print(ggplot(data = PAC_list[[PAC]][["high_load"]], aes(x = omnibus_span_no_DFR, y = thetP_lgamA_mnt))+
          geom_point()+
          stat_smooth(method="lm")+
          ggtitle(names(PAC_list)[PAC]))
  
  print(ggplot(data = PAC_list[[PAC]][["high_load"]], aes(x = omnibus_span_no_DFR, y = thetP_hgamA_mnt))+
          geom_point()+
          stat_smooth(method="lm")+
          ggtitle(names(PAC_list)[PAC]))
  
  print(ggplot(data = PAC_list[[PAC]][["high_load"]], aes(x = omnibus_span_no_DFR, y = thetP_ugamA_mnt))+
          geom_point()+
          stat_smooth(method="lm")+
          ggtitle(names(PAC_list)[PAC]))
  
  
}

Load Effect

for (PAC in seq.int(1,length(PAC_list))){
  
  print(ggplot(data = PAC_list[[PAC]][["load_effect"]], aes(x = omnibus_span_no_DFR, y = thetP_hbetA_mnt))+
          geom_point()+
          stat_smooth(method="lm")+
          ggtitle(names(PAC_list)[PAC]))
  
  print(ggplot(data = PAC_list[[PAC]][["load_effect"]], aes(x = omnibus_span_no_DFR, y = thetP_lgamA_mnt))+
          geom_point()+
          stat_smooth(method="lm")+
          ggtitle(names(PAC_list)[PAC]))
  
  print(ggplot(data = PAC_list[[PAC]][["load_effect"]], aes(x = omnibus_span_no_DFR, y = thetP_hgamA_mnt))+
          geom_point()+
          stat_smooth(method="lm")+
          ggtitle(names(PAC_list)[PAC]))
  
  print(ggplot(data = PAC_list[[PAC]][["load_effect"]], aes(x = omnibus_span_no_DFR, y = thetP_ugamA_mnt))+
          geom_point()+
          stat_smooth(method="lm")+
          ggtitle(names(PAC_list)[PAC]))
  
  
}

Alpha - Gamma

High Load

for (PAC in seq.int(1,length(PAC_list))){
  
  print(ggplot(data = PAC_list[[PAC]][["high_load"]], aes(x = omnibus_span_no_DFR, y = alphP_lgamA_mnt))+
          geom_point()+
          stat_smooth(method="lm")+
          ggtitle(names(PAC_list)[PAC]))
  
  print(ggplot(data = PAC_list[[PAC]][["high_load"]], aes(x = omnibus_span_no_DFR, y = alphP_hgamA_mnt))+
          geom_point()+
          stat_smooth(method="lm")+
          ggtitle(names(PAC_list)[PAC]))
  
  print(ggplot(data = PAC_list[[PAC]][["high_load"]], aes(x = omnibus_span_no_DFR, y = alphP_ugamA_mnt))+
          geom_point()+
          stat_smooth(method="lm")+
          ggtitle(names(PAC_list)[PAC]))
  
  
}

## Warning: Removed 17 rows containing non-finite values (stat_smooth).

## Warning: Removed 17 rows containing missing values (geom_point).

## Warning: Removed 17 rows containing non-finite values (stat_smooth).

## Warning: Removed 17 rows containing missing values (geom_point).

## Warning: Removed 17 rows containing non-finite values (stat_smooth).

## Warning: Removed 17 rows containing missing values (geom_point).

## Warning: Removed 27 rows containing non-finite values (stat_smooth).

## Warning: Removed 27 rows containing missing values (geom_point).

## Warning: Removed 27 rows containing non-finite values (stat_smooth).

## Warning: Removed 27 rows containing missing values (geom_point).

## Warning: Removed 27 rows containing non-finite values (stat_smooth).

## Warning: Removed 27 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

Load Effects

for (PAC in seq.int(1,length(PAC_list))){
  
  print(ggplot(data = PAC_list[[PAC]][["load_effect"]], aes(x = omnibus_span_no_DFR, y = alphP_lgamA_mnt))+
          geom_point()+
          stat_smooth(method="lm")+
          ggtitle(names(PAC_list)[PAC]))
  
  print(ggplot(data = PAC_list[[PAC]][["load_effect"]], aes(x = omnibus_span_no_DFR, y = alphP_hgamA_mnt))+
          geom_point()+
          stat_smooth(method="lm")+
          ggtitle(names(PAC_list)[PAC]))
  
  print(ggplot(data = PAC_list[[PAC]][["load_effect"]], aes(x = omnibus_span_no_DFR, y = alphP_ugamA_mnt))+
          geom_point()+
          stat_smooth(method="lm")+
          ggtitle(names(PAC_list)[PAC]))
  
  
}

## Warning: Removed 17 rows containing non-finite values (stat_smooth).

## Warning: Removed 17 rows containing missing values (geom_point).

## Warning: Removed 17 rows containing non-finite values (stat_smooth).

## Warning: Removed 17 rows containing missing values (geom_point).

## Warning: Removed 17 rows containing non-finite values (stat_smooth).

## Warning: Removed 17 rows containing missing values (geom_point).

## Warning: Removed 27 rows containing non-finite values (stat_smooth).

## Warning: Removed 27 rows containing missing values (geom_point).

## Warning: Removed 27 rows containing non-finite values (stat_smooth).

## Warning: Removed 27 rows containing missing values (geom_point).

## Warning: Removed 27 rows containing non-finite values (stat_smooth).

## Warning: Removed 27 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

## Warning: Removed 7 rows containing non-finite values (stat_smooth).

## Warning: Removed 7 rows containing missing values (geom_point).

Bar plots

split_PAC <- list()
split_PAC_prepped <- list()

split_PAC_high <- list()
split_PAC_high_prepped <- list()

for (PAC in seq.int(1,length(PAC_list))){
  split_PAC[[PAC]] <- split_into_groups(PAC_list[[PAC]][["load_effect"]], WM_groups)
  split_PAC_prepped[[PAC]] <- prep_split_for_bar_plots(split_PAC[[PAC]])
  
  split_PAC_high[[PAC]] <- split_into_groups(PAC_list[[PAC]][["high_load"]], WM_groups)
  split_PAC_high_prepped[[PAC]] <- prep_split_for_bar_plots(split_PAC_high[[PAC]])
}

names(split_PAC) <- names(PAC_list)
names(split_PAC_prepped) <- names(PAC_list)

names(split_PAC_high) <- names(PAC_list)
names(split_PAC_high_prepped) <- names(PAC_list)

PAC_bar_plots_high <- list(list(), list(),list(), list(),list(), list())

for (PAC in seq.int(1,length(PAC_list))){
  split_PAC_high_prepped[[PAC]][["melt_data"]] %>% 
    filter(variable == "thetP_hbetA_mnt") %>% 
    mutate(level = factor(level, levels =c("low","med","high"))) %>%
    ggplot()+
    geom_bar(aes(x=level,y=Means),stat="identity")+
    geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
    #ggtitle(names(PAC_list)[PAC])+ 
    xlab("WMC")+
    ylab("thetP_betA_mnt")+
    theme_classic()+
    theme(aspect.ratio = 1) -> PAC_bar_plots_high[[PAC]][["thetP_betA_mnt"]]
  
  split_PAC_high_prepped[[PAC]][["melt_data"]] %>% 
    filter(variable == "thetP_lgamA_mnt") %>% 
    mutate(level = factor(level, levels =c("low","med","high"))) %>%
    ggplot()+
    geom_bar(aes(x=level,y=Means),stat="identity")+
    geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
    #ggtitle(names(PAC_list)[PAC])+ 
    xlab("WMC")+
    ylab("thetP_lgamA_mnt")+
    theme_classic()+
    theme(aspect.ratio = 1) -> PAC_bar_plots_high[[PAC]][["thetP_lgamA_mnt"]]
  
  split_PAC_high_prepped[[PAC]][["melt_data"]] %>% 
    filter(variable == "thetP_hgamA_mnt") %>% 
    mutate(level = factor(level, levels =c("low","med","high"))) %>%
    ggplot()+
    geom_bar(aes(x=level,y=Means),stat="identity")+
    geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
    #ggtitle(names(PAC_list)[PAC])+ 
    xlab("WMC")+
    ylab("thetP_hgamA_mnt")+
    theme_classic()+
    theme(aspect.ratio = 1) -> PAC_bar_plots_high[[PAC]][["thetP_hgamA_mnt"]]
  
  split_PAC_high_prepped[[PAC]][["melt_data"]] %>% 
    filter(variable == "thetP_ugamA_mnt") %>% 
    mutate(level = factor(level, levels =c("low","med","high"))) %>%
    ggplot()+
    geom_bar(aes(x=level,y=Means),stat="identity")+
    geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
    #ggtitle(names(PAC_list)[PAC])+ 
    xlab("WMC")+
    ylab("thetP_ugamA_mnt")+
    theme_classic()+
    theme(aspect.ratio = 1) -> PAC_bar_plots_high[[PAC]][["thetP_ugamA_mnt"]]
  
}

names(PAC_bar_plots_high) <- names(PAC_list)

PAC_bar_plots <- list(list(), list(),list(), list(),list(), list())

for (PAC in seq.int(1,length(PAC_list))){
  split_PAC_prepped[[PAC]][["melt_data"]] %>% 
    filter(variable == "thetP_hbetA_mnt") %>% 
    mutate(level = factor(level, levels =c("low","med","high"))) %>%
    ggplot()+
    geom_bar(aes(x=level,y=Means),stat="identity")+
    geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
    #ggtitle(names(PAC_list)[PAC])+ 
    xlab("WMC")+
    ylab("thetP_betA_mnt")+
    theme_classic()+
    theme(aspect.ratio = 1) -> PAC_bar_plots[[PAC]][["thetP_betA_mnt"]]
  
  split_PAC_prepped[[PAC]][["melt_data"]] %>% 
    filter(variable == "thetP_lgamA_mnt") %>% 
    mutate(level = factor(level, levels =c("low","med","high"))) %>%
    ggplot()+
    geom_bar(aes(x=level,y=Means),stat="identity")+
    geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
    #ggtitle(names(PAC_list)[PAC])+ 
    xlab("WMC")+
    ylab("thetP_lgamA_mnt")+
    theme_classic()+
    theme(aspect.ratio = 1) -> PAC_bar_plots[[PAC]][["thetP_lgamA_mnt"]]
  
  split_PAC_prepped[[PAC]][["melt_data"]] %>% 
    filter(variable == "thetP_hgamA_mnt") %>% 
    mutate(level = factor(level, levels =c("low","med","high"))) %>%
    ggplot()+
    geom_bar(aes(x=level,y=Means),stat="identity")+
    geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
    #ggtitle(names(PAC_list)[PAC])+ 
    xlab("WMC")+
    ylab("thetP_hgamA_mnt")+
    theme_classic()+
    theme(aspect.ratio = 1) -> PAC_bar_plots[[PAC]][["thetP_hgamA_mnt"]]
  
  split_PAC_prepped[[PAC]][["melt_data"]] %>% 
    filter(variable == "thetP_ugamA_mnt") %>% 
    mutate(level = factor(level, levels =c("low","med","high"))) %>%
    ggplot()+
    geom_bar(aes(x=level,y=Means),stat="identity")+
    geom_errorbar(aes(x=level,ymin=Means-SE, ymax=Means+SE), width=0.2)+
    #ggtitle(names(PAC_list)[PAC])+ 
    xlab("WMC")+
    ylab("thetP_ugamA_mnt")+
    theme_classic()+
    theme(aspect.ratio = 1) -> PAC_bar_plots[[PAC]][["thetP_ugamA_mnt"]]
  
}

names(PAC_bar_plots) <- names(PAC_list)

for (PAC in seq.int(1,length(PAC_bar_plots_high))){
  
  print(
    (PAC_bar_plots_high[[PAC]][[1]] + PAC_bar_plots_high[[PAC]][[2]]) / 
      (PAC_bar_plots_high[[PAC]][[3]] + PAC_bar_plots_high[[PAC]][[4]]) +
      plot_annotation(title = names(PAC_bar_plots_high)[PAC])
  )
  
}

for (PAC in seq.int(1,length(PAC_list))){
  for (measure in seq.int(2,5)){
    print(paste0("Region: ", names(PAC_list)[PAC], ", measure: ",colnames(PAC_list[[PAC]][["high_load"]])[measure]))
    temp.aov <-  aov(PAC_list[[PAC]][["high_load"]][,measure] ~ PAC_list[[PAC]][["high_load"]]$level)
    print(summary(temp.aov))
  }
}

## [1] "Region: midOccip, measure: thetP_hbetA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   2.35   1.174   0.717   0.49
## Residuals                            148 242.17   1.636               
## 17 observations deleted due to missingness
## [1] "Region: midOccip, measure: thetP_lgamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   2.17   1.087   0.878  0.418
## Residuals                            148 183.38   1.239               
## 17 observations deleted due to missingness
## [1] "Region: midOccip, measure: thetP_hgamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   0.62  0.3087   1.298  0.276
## Residuals                            148  35.21  0.2379               
## 17 observations deleted due to missingness
## [1] "Region: midOccip, measure: thetP_ugamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   0.55  0.2748   1.317  0.271
## Residuals                            148  30.89  0.2087               
## 17 observations deleted due to missingness
## [1] "Region: midFront, measure: thetP_hbetA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   0.81  0.4053   1.125  0.328
## Residuals                            138  49.72  0.3603               
## 27 observations deleted due to missingness
## [1] "Region: midFront, measure: thetP_lgamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   0.15 0.07282   0.245  0.783
## Residuals                            138  40.93 0.29660               
## 27 observations deleted due to missingness
## [1] "Region: midFront, measure: thetP_hgamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   0.08 0.04154   0.163   0.85
## Residuals                            138  35.13 0.25455               
## 27 observations deleted due to missingness
## [1] "Region: midFront, measure: thetP_ugamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   0.19 0.09441   0.385  0.681
## Residuals                            138  33.86 0.24539               
## 27 observations deleted due to missingness
## [1] "Region: Oz, measure: thetP_hbetA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   12.6   6.322   1.862  0.159
## Residuals                            158  536.6   3.396               
## 7 observations deleted due to missingness
## [1] "Region: Oz, measure: thetP_lgamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   0.19  0.0946   0.074  0.929
## Residuals                            158 202.42  1.2812               
## 7 observations deleted due to missingness
## [1] "Region: Oz, measure: thetP_hgamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   0.09  0.0442   0.044  0.957
## Residuals                            158 159.59  1.0101               
## 7 observations deleted due to missingness
## [1] "Region: Oz, measure: thetP_ugamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   1.71  0.8554   0.573  0.565
## Residuals                            158 235.89  1.4930               
## 7 observations deleted due to missingness
## [1] "Region: O1, measure: thetP_hbetA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2    0.3  0.1618    0.06  0.942
## Residuals                            158  427.1  2.7033               
## 7 observations deleted due to missingness
## [1] "Region: O1, measure: thetP_lgamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2    1.8  0.9228   0.414  0.662
## Residuals                            158  352.1  2.2286               
## 7 observations deleted due to missingness
## [1] "Region: O1, measure: thetP_hgamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   0.37   0.187   0.184  0.832
## Residuals                            158 160.90   1.018               
## 7 observations deleted due to missingness
## [1] "Region: O1, measure: thetP_ugamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   3.36   1.678   1.384  0.254
## Residuals                            158 191.51   1.212               
## 7 observations deleted due to missingness
## [1] "Region: O2, measure: thetP_hbetA_mnt"
##                                       Df Sum Sq Mean Sq F value  Pr(>F)   
## PAC_list[[PAC]][["high_load"]]$level   2   29.5  14.752   5.016 0.00772 **
## Residuals                            158  464.7   2.941                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 7 observations deleted due to missingness
## [1] "Region: O2, measure: thetP_lgamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2    1.2  0.6008    0.49  0.614
## Residuals                            158  193.9  1.2270               
## 7 observations deleted due to missingness
## [1] "Region: O2, measure: thetP_hgamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2    1.0  0.5025   0.382  0.683
## Residuals                            158  208.1  1.3170               
## 7 observations deleted due to missingness
## [1] "Region: O2, measure: thetP_ugamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   8.83   4.417   2.208  0.113
## Residuals                            158 315.99   2.000               
## 7 observations deleted due to missingness
## [1] "Region: Fz, measure: thetP_hbetA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)  
## PAC_list[[PAC]][["high_load"]]$level   2   3.43  1.7170   2.916 0.0571 .
## Residuals                            158  93.03  0.5888                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 7 observations deleted due to missingness
## [1] "Region: Fz, measure: thetP_lgamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   1.38  0.6886   0.927  0.398
## Residuals                            158 117.40  0.7430               
## 7 observations deleted due to missingness
## [1] "Region: Fz, measure: thetP_hgamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)  
## PAC_list[[PAC]][["high_load"]]$level   2   1.92  0.9615   2.959 0.0548 .
## Residuals                            158  51.34  0.3250                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 7 observations deleted due to missingness
## [1] "Region: Fz, measure: thetP_ugamA_mnt"
##                                       Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["high_load"]]$level   2   0.53  0.2651   0.707  0.495
## Residuals                            158  59.27  0.3751               
## 7 observations deleted due to missingness

for (PAC in seq.int(1,length(PAC_bar_plots))){
  
  print(
    (PAC_bar_plots[[PAC]][[1]] + PAC_bar_plots[[PAC]][[2]]) / 
      (PAC_bar_plots[[PAC]][[3]] + PAC_bar_plots[[PAC]][[4]]) +
      plot_annotation(title = names(PAC_bar_plots)[PAC])
  )
  
}

for (PAC in seq.int(1,length(PAC_list))){
  for (measure in seq.int(2,5)){
    print(paste0("Region: ", names(PAC_list)[PAC], ", measure: ",colnames(PAC_list[[PAC]][["load_effect"]])[measure]))
    temp.aov <-  aov(PAC_list[[PAC]][["load_effect"]][,measure] ~ PAC_list[[PAC]][["load_effect"]]$level)
    print(summary(temp.aov))
  }
}

## [1] "Region: midOccip, measure: thetP_hbetA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   5.46   2.731   1.952  0.146
## Residuals                              148 207.02   1.399               
## 17 observations deleted due to missingness
## [1] "Region: midOccip, measure: thetP_lgamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   3.44   1.721    1.52  0.222
## Residuals                              148 167.57   1.132               
## 17 observations deleted due to missingness
## [1] "Region: midOccip, measure: thetP_hgamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   0.90   0.452   0.855  0.428
## Residuals                              148  78.29   0.529               
## 17 observations deleted due to missingness
## [1] "Region: midOccip, measure: thetP_ugamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   0.06  0.0297   0.061  0.941
## Residuals                              148  71.94  0.4861               
## 17 observations deleted due to missingness
## [1] "Region: midFront, measure: thetP_hbetA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   1.22  0.6088   1.076  0.344
## Residuals                              138  78.05  0.5656               
## 27 observations deleted due to missingness
## [1] "Region: midFront, measure: thetP_lgamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   1.05  0.5239   1.001   0.37
## Residuals                              138  72.25  0.5235               
## 27 observations deleted due to missingness
## [1] "Region: midFront, measure: thetP_hgamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   0.18  0.0890   0.194  0.824
## Residuals                              138  63.41  0.4595               
## 27 observations deleted due to missingness
## [1] "Region: midFront, measure: thetP_ugamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   0.12  0.0586   0.119  0.888
## Residuals                              138  67.80  0.4913               
## 27 observations deleted due to missingness
## [1] "Region: Oz, measure: thetP_hbetA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2    9.0   4.490   1.349  0.263
## Residuals                              158  525.9   3.329               
## 7 observations deleted due to missingness
## [1] "Region: Oz, measure: thetP_lgamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   2.84   1.421   0.753  0.473
## Residuals                              158 298.31   1.888               
## 7 observations deleted due to missingness
## [1] "Region: Oz, measure: thetP_hgamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   4.36   2.178   1.159  0.316
## Residuals                              158 296.85   1.879               
## 7 observations deleted due to missingness
## [1] "Region: Oz, measure: thetP_ugamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   6.23   3.115   1.748  0.177
## Residuals                              158 281.57   1.782               
## 7 observations deleted due to missingness
## [1] "Region: O1, measure: thetP_hbetA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2    4.1   2.064   0.838  0.435
## Residuals                              158  389.2   2.463               
## 7 observations deleted due to missingness
## [1] "Region: O1, measure: thetP_lgamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2    0.4  0.2044   0.087  0.917
## Residuals                              158  371.2  2.3495               
## 7 observations deleted due to missingness
## [1] "Region: O1, measure: thetP_hgamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   0.12  0.0622   0.045  0.956
## Residuals                              158 219.13  1.3869               
## 7 observations deleted due to missingness
## [1] "Region: O1, measure: thetP_ugamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2    0.8   0.397   0.178  0.837
## Residuals                              158  353.1   2.235               
## 7 observations deleted due to missingness
## [1] "Region: O2, measure: thetP_hbetA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)  
## PAC_list[[PAC]][["load_effect"]]$level   2   12.2   6.103   2.396 0.0944 .
## Residuals                              158  402.5   2.548                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 7 observations deleted due to missingness
## [1] "Region: O2, measure: thetP_lgamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2    4.4   2.206   0.956  0.387
## Residuals                              158  364.5   2.307               
## 7 observations deleted due to missingness
## [1] "Region: O2, measure: thetP_hgamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2    2.1   1.071   0.358    0.7
## Residuals                              158  473.1   2.994               
## 7 observations deleted due to missingness
## [1] "Region: O2, measure: thetP_ugamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)  
## PAC_list[[PAC]][["load_effect"]]$level   2   13.1   6.572   2.505 0.0849 .
## Residuals                              158  414.4   2.623                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 7 observations deleted due to missingness
## [1] "Region: Fz, measure: thetP_hbetA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)  
## PAC_list[[PAC]][["load_effect"]]$level   2   9.58   4.790   4.522 0.0123 *
## Residuals                              158 167.36   1.059                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 7 observations deleted due to missingness
## [1] "Region: Fz, measure: thetP_lgamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   2.03   1.016   0.987  0.375
## Residuals                              158 162.73   1.030               
## 7 observations deleted due to missingness
## [1] "Region: Fz, measure: thetP_hgamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   0.47  0.2359   0.373   0.69
## Residuals                              158 100.03  0.6331               
## 7 observations deleted due to missingness
## [1] "Region: Fz, measure: thetP_ugamA_mnt"
##                                         Df Sum Sq Mean Sq F value Pr(>F)
## PAC_list[[PAC]][["load_effect"]]$level   2   1.04  0.5225   0.715  0.491
## Residuals                              158 115.46  0.7308               
## 7 observations deleted due to missingness

CDA

CDA <- data.frame(read.csv("data/EEG/LCD/ERPS_CDA.txt"))
CDA_fMRI <- CDA[CDA$subID %in% constructs_fMRI$PTID,]
CDA_fMRI[is.nan(CDA_fMRI[,2]),2:11] <- NA
colnames(CDA_fMRI)[1] <- "PTID"
CDA_fMRI$CDA_3_1 <- CDA_fMRI$LCD_CDA_3dots - CDA_fMRI$LCD_CDA_1dots
CDA_fMRI$CDA_5_1 <- CDA_fMRI$LCD_CDA_5dots - CDA_fMRI$LCD_CDA_1dots
CDA_fMRI$CDA_5_3 <- CDA_fMRI$LCD_CDA_5dots - CDA_fMRI$LCD_CDA_3dots

Raw Data

First, let’s take a look at the raw data, just for sanity’s sake.

CDA_fMRI %>% 
  select(c("LCD_CDA_1dots", "LCD_CDA_3dots", "LCD_CDA_5dots")) %>% 
  melt() %>% 
  ggplot(aes(x=variable,y=value))+
    geom_jitter()+
    stat_summary(fun.data=mean_sdl, geom="pointrange", color="red")+
    xlab("Level")+
  ylab("CDA")

## No id variables; using all as measure variables

## Warning: Removed 24 rows containing non-finite values (stat_summary).

## Warning: Removed 24 rows containing missing values (geom_point).

CDA_fMRI %>% 
  select(c("CDA_3_1", "CDA_5_1", "CDA_5_3")) %>% 
  melt() %>% 
  ggplot(aes(x=variable,y=value))+
    geom_jitter()+
    stat_summary(fun.data=mean_sdl, geom="pointrange", color="red")+
    xlab("Level differences")+
  ylab("CDA Load Effects")

## No id variables; using all as measure variables

## Warning: Removed 24 rows containing non-finite values (stat_summary).

## Warning: Removed 24 rows containing missing values (geom_point).

Relationship of Raw Data with Span

Other studies have shown a relationship between the CDA and capacity as measured by LCD K. Our data, only showed a correlation between the CDA in the low load condition (1 dot) and the LCD K max, with no correlation between omnibus span or DFR load effect.

data_to_plot <- data.frame(CDA_fMRI, p200_delay_DFR, omnibus_span = constructs_fMRI$omnibus_span_no_DFR)
data_to_plot <- merge(data_to_plot,p200_cog)

(ggplot(data_to_plot,aes(x=DFR_Load3_Load1, y=LCD_CDA_5dots))+
    geom_point()+
    stat_smooth(method="loess")+
    theme_classic())+(ggplot(data_to_plot,aes(x=DFR_Load3_Load1, y=LCD_CDA_3dots))+
                        geom_point()+
                        stat_smooth(method="lm")+
                        theme_classic())+(ggplot(data_to_plot,aes(x=DFR_Load3_Load1, y=LCD_CDA_1dots))+
                                            geom_point()+
                                            stat_smooth(method="lm")+
                                            theme_classic()+
                                            theme(aspect.ratio=1))

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

cor.test(data_to_plot$LCD_CDA_5dots, data_to_plot$DFR_Load3_Load1)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_5dots and data_to_plot$DFR_Load3_Load1
## t = 0.1352, df = 160, p-value = 0.8926
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1437443  0.1646122
## sample estimates:
##        cor 
## 0.01068802

cor.test(data_to_plot$LCD_CDA_3dots, data_to_plot$DFR_Load3_Load1)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_3dots and data_to_plot$DFR_Load3_Load1
## t = -1.524, df = 160, p-value = 0.1295
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.26885513  0.03522659
## sample estimates:
##        cor 
## -0.1196186

cor.test(data_to_plot$LCD_CDA_1dots, data_to_plot$DFR_Load3_Load1)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_1dots and data_to_plot$DFR_Load3_Load1
## t = 0.92497, df = 160, p-value = 0.3564
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.08218942  0.22460001
## sample estimates:
##       cor 
## 0.0729303

(ggplot(data_to_plot,aes(x=omnibus_span, y=LCD_CDA_5dots))+
    geom_point()+
    stat_smooth(method="loess")+
    theme_classic())+(ggplot(data_to_plot,aes(x=omnibus_span, y=LCD_CDA_3dots))+
                        geom_point()+
                        stat_smooth(method="lm")+
                        theme_classic())+(ggplot(data_to_plot,aes(x=omnibus_span, y=LCD_CDA_1dots))+
                                            geom_point()+
                                            stat_smooth(method="lm")+
                                            theme_classic()+
                                            theme(aspect.ratio=1))

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

cor.test(data_to_plot$LCD_CDA_5dots, data_to_plot$omnibus_span)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_5dots and data_to_plot$omnibus_span
## t = -0.81334, df = 160, p-value = 0.4172
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.21622386  0.09092732
## sample estimates:
##         cor 
## -0.06416781

cor.test(data_to_plot$LCD_CDA_3dots, data_to_plot$omnibus_span)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_3dots and data_to_plot$omnibus_span
## t = -1.4802, df = 160, p-value = 0.1408
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.26566128  0.03866141
## sample estimates:
##        cor 
## -0.1162269

cor.test(data_to_plot$LCD_CDA_1dots, data_to_plot$omnibus_span)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_1dots and data_to_plot$omnibus_span
## t = -1.3219, df = 160, p-value = 0.1881
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.25406112  0.05107662
## sample estimates:
##        cor 
## -0.1039374

(ggplot(data_to_plot,aes(x=XLCD_K_MAX, y=LCD_CDA_5dots))+
    geom_point()+
    stat_smooth(method="loess")+
    theme_classic())+(ggplot(data_to_plot,aes(x=XLCD_K_MAX, y=LCD_CDA_3dots))+
                        geom_point()+
                        stat_smooth(method="lm")+
                        theme_classic())+(ggplot(data_to_plot,aes(x=XLCD_K_MAX, y=LCD_CDA_1dots))+
                                            geom_point()+
                                            stat_smooth(method="lm")+
                                            theme_classic()+
                                            theme(aspect.ratio=1))

## Warning: Removed 9 rows containing non-finite values (stat_smooth).

## Warning: Removed 9 rows containing missing values (geom_point).

## Warning: Removed 9 rows containing non-finite values (stat_smooth).

## Warning: Removed 9 rows containing missing values (geom_point).

## Warning: Removed 9 rows containing non-finite values (stat_smooth).

## Warning: Removed 9 rows containing missing values (geom_point).

cor.test(data_to_plot$LCD_CDA_5dots, data_to_plot$XLCD_K_MAX)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_5dots and data_to_plot$XLCD_K_MAX
## t = -0.48689, df = 159, p-value = 0.627
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1921126  0.1167875
## sample estimates:
##         cor 
## -0.03858431

cor.test(data_to_plot$LCD_CDA_3dots, data_to_plot$XLCD_K_MAX)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_3dots and data_to_plot$XLCD_K_MAX
## t = -1.2834, df = 159, p-value = 0.2012
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.25198309  0.05427012
## sample estimates:
##        cor 
## -0.1012547

cor.test(data_to_plot$LCD_CDA_1dots, data_to_plot$XLCD_K_MAX)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$LCD_CDA_1dots and data_to_plot$XLCD_K_MAX
## t = -2.0777, df = 159, p-value = 0.03935
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.309470392 -0.008106698
## sample estimates:
##        cor 
## -0.1625777

CDA_split <- split_into_groups(CDA_fMRI, WM_groups)
CDA_split_plot_data <- prep_split_for_bar_plots(CDA_split)

CDA_split_plot_data[["melt_data"]] %>%
  filter(variable == "LCD_CDA_5dots") %>%
  mutate(level = factor(level,levels=c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE,ymax=Means+SE), width=0.2)+
  xlab("WMC group")+
  ylab("CDA")+
  ggtitle("CDA - 5 dots")+ 
  theme_classic()+
  theme(aspect.ratio=1)

CDA_split_plot_data[["melt_data"]] %>%
  filter(variable == "LCD_CDA_1dots") %>%
  mutate(level = factor(level,levels=c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE,ymax=Means+SE), width=0.2)+
  xlab("WMC group")+
  ylab("CDA")+
  ggtitle("CDA - 1 dot")+ 
  theme_classic()+
  theme(aspect.ratio=1)

CDA_anova <- merge(CDA_fMRI,WM_to_merge, by="PTID")

summary(aov(LCD_CDA_1dots ~ level,data=CDA_anova))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.25  0.1266    0.43  0.651
## Residuals   157  46.25  0.2946               
## 8 observations deleted due to missingness

Relationship of Load Effects with Span

Vogel and Machizawa (2004) had mean K = 2.8. Our average K max = 2.68 (so slightly lower). The only thing that even shows a trend with DFR fMRI load effects is in the load effect between the 3 dot and 1 dot condition.

(ggplot(data_to_plot,aes(x=DFR_Load3_Load1, y=CDA_5_3))+
   geom_point()+
   stat_smooth(method="loess")+
   theme_classic())+(ggplot(data_to_plot,aes(x=DFR_Load3_Load1, y=CDA_5_1))+
                       geom_point()+
                       stat_smooth(method="lm")+
                       theme_classic())+(ggplot(data_to_plot,aes(x=DFR_Load3_Load1, y=CDA_3_1))+
                                           geom_point()+
                                           stat_smooth(method="lm")+
                                           theme_classic()+
                                           theme(aspect.ratio=1))

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

cor.test(data_to_plot$CDA_5_3, data_to_plot$DFR_Load3_Load1)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_5_3 and data_to_plot$DFR_Load3_Load1
## t = 1.3562, df = 160, p-value = 0.177
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.04838862  0.25658013
## sample estimates:
##       cor 
## 0.1066022

cor.test(data_to_plot$CDA_5_1, data_to_plot$DFR_Load3_Load1)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_5_1 and data_to_plot$DFR_Load3_Load1
## t = -0.56408, df = 160, p-value = 0.5735
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1973896  0.1104038
## sample estimates:
##         cor 
## -0.04455009

cor.test(data_to_plot$CDA_3_1, data_to_plot$DFR_Load3_Load1)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_3_1 and data_to_plot$DFR_Load3_Load1
## t = -1.9681, df = 160, p-value = 0.05079
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3008041900  0.0004669451
## sample estimates:
##        cor 
## -0.1537396

(ggplot(data_to_plot,aes(x=omnibus_span, y=CDA_5_3))+
    geom_point()+
    stat_smooth(method="lm")+
    theme_classic())+(ggplot(data_to_plot,aes(x=omnibus_span, y=CDA_5_1))+
                        geom_point()+
                        stat_smooth(method="lm")+
                        theme_classic())+(ggplot(data_to_plot,aes(x=omnibus_span, y=CDA_3_1))+
                                            geom_point()+
                                            stat_smooth(method="lm")+
                                            theme_classic()+
                                            theme(aspect.ratio=1)
                        )

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

## Warning: Removed 8 rows containing non-finite values (stat_smooth).

## Warning: Removed 8 rows containing missing values (geom_point).

cor.test(data_to_plot$CDA_5_3, data_to_plot$omnibus_span)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_5_3 and data_to_plot$omnibus_span
## t = 0.59677, df = 160, p-value = 0.5515
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1078529  0.1998694
## sample estimates:
##        cor 
## 0.04712628

cor.test(data_to_plot$CDA_5_1, data_to_plot$omnibus_span)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_5_1 and data_to_plot$omnibus_span
## t = 0.30731, df = 160, p-value = 0.759
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1303964  0.1778170
## sample estimates:
##        cor 
## 0.02428743

cor.test(data_to_plot$CDA_3_1, data_to_plot$omnibus_span)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_3_1 and data_to_plot$omnibus_span
## t = -0.31616, df = 160, p-value = 0.7523
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1784947  0.1297082
## sample estimates:
##       cor 
## -0.024987

mean(data_to_plot$XLCD_K_MAX,na.rm=TRUE)

## [1] 2.686799

(ggplot(data_to_plot,aes(x=XLCD_K_MAX, y=CDA_5_3))+
    geom_point()+
    stat_smooth(method="loess")+
    theme_classic()+
    theme(aspect.ratio=1))+(ggplot(data_to_plot,aes(x=XLCD_K_MAX, y=CDA_5_1))+
                              geom_point()+
                              stat_smooth(method="loess")+
                              theme_classic())+(ggplot(data_to_plot,aes(x=XLCD_K_MAX, y=CDA_3_1))+
                                                  geom_point()+
                                                  stat_smooth(method="loess")+
                                                  theme_classic()) 

## Warning: Removed 9 rows containing non-finite values (stat_smooth).

## Warning: Removed 9 rows containing missing values (geom_point).

## Warning: Removed 9 rows containing non-finite values (stat_smooth).

## Warning: Removed 9 rows containing missing values (geom_point).

## Warning: Removed 9 rows containing non-finite values (stat_smooth).

## Warning: Removed 9 rows containing missing values (geom_point).

cor.test(data_to_plot$CDA_5_3, data_to_plot$XLCD_K_MAX)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_5_3 and data_to_plot$XLCD_K_MAX
## t = 0.68698, df = 159, p-value = 0.4931
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1011252  0.2073309
## sample estimates:
##        cor 
## 0.05440059

cor.test(data_to_plot$CDA_5_1, data_to_plot$XLCD_K_MAX)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_5_1 and data_to_plot$XLCD_K_MAX
## t = 1.118, df = 159, p-value = 0.2652
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.06727587  0.23971818
## sample estimates:
##        cor 
## 0.08831801

cor.test(data_to_plot$CDA_3_1, data_to_plot$XLCD_K_MAX)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$CDA_3_1 and data_to_plot$XLCD_K_MAX
## t = 0.37401, df = 159, p-value = 0.7089
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.1256028  0.1834815
## sample estimates:
##        cor 
## 0.02964807

ggplot(data_to_plot,aes(x=XLCD_K_MAX, y=omnibus_span))+
  geom_point()+
  stat_smooth(method="lm")+
  theme_classic()

## Warning: Removed 3 rows containing non-finite values (stat_smooth).

## Warning: Removed 3 rows containing missing values (geom_point).

cor.test(data_to_plot$omnibus_span, data_to_plot$XLCD_K_MAX)

## 
##  Pearson's product-moment correlation
## 
## data:  data_to_plot$omnibus_span and data_to_plot$XLCD_K_MAX
## t = 8.8404, df = 165, p-value = 1.384e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.4541707 0.6618156
## sample estimates:
##       cor 
## 0.5669319

Numerically, it looks like there is an inverted U shaped relationship for the 5 dots - 1 dot and 3 dot - 1 dot, but those relationships are not significant.

CDA_split_plot_data[["melt_data"]] %>%
  filter(variable == "CDA_5_3") %>%
  mutate(level = factor(level,levels=c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE,ymax=Means+SE), width=0.2)+
  xlab("WMC group")+
  ylab("CDA Load Effect")+
  ggtitle("5 dots - 3 dots")+ 
  theme_classic()+
  theme(aspect.ratio=1)

CDA_split_plot_data[["melt_data"]] %>%
  filter(variable == "CDA_5_1") %>%
  mutate(level = factor(level,levels=c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE,ymax=Means+SE), width=0.2)+
  xlab("WMC group")+
  ylab("CDA Load Effect")+
  ggtitle("5 dots - 1 dot")+ 
  theme_classic()+
  theme(aspect.ratio=1)

CDA_split_plot_data[["melt_data"]] %>%
  filter(variable == "CDA_3_1") %>%
  mutate(level = factor(level,levels=c("low","med","high"))) %>%
  ggplot()+
  geom_bar(aes(x=level,y=Means),stat="identity")+
  geom_errorbar(aes(x=level,ymin=Means-SE,ymax=Means+SE), width=0.2)+
  xlab("WMC group")+
  ylab("CDA Load Effect")+
  ggtitle("3 dots - 1 dot")+ 
  theme_classic()+
  theme(aspect.ratio=1)

CDA_5_3.aov <- aov(CDA_5_3 ~ level, data = CDA_anova)
print("CDA - 5 dots - 3 dots")

## [1] "CDA - 5 dots - 3 dots"

print(summary(CDA_5_3.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2    1.0  0.4984   0.782  0.459
## Residuals   157  100.1  0.6376               
## 8 observations deleted due to missingness

CDA_5_1.aov <- aov(CDA_5_1 ~ level, data = CDA_anova)
print("CDA - 5 dots - 1 dot")

## [1] "CDA - 5 dots - 1 dot"

print(summary(CDA_5_1.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   0.33  0.1675     0.3  0.741
## Residuals   157  87.73  0.5588               
## 8 observations deleted due to missingness

CDA_3_1.aov <- aov(CDA_3_1 ~ level, data = CDA_anova)
print("CDA - 3 dots - 1 dot")

## [1] "CDA - 3 dots - 1 dot"

print(summary(CDA_3_1.aov))

##              Df Sum Sq Mean Sq F value Pr(>F)
## level         2   1.13  0.5650   0.958  0.386
## Residuals   157  92.59  0.5897               
## 8 observations deleted due to missingness